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The  objective  of  this  research  is  twofold.  One  is  to  develop  a parametric  model  of  the 
pulse-wave  ultrasound  Doppler  signal.  The  other  is  to  propose  a blood  velocity 
estimation  technique  based  on  the  autoregressive  spectral  estimation  (ARSE). 

The  spectral  analysis  of  the  Doppler  signal  provides  the  velocity  information  in  the 
sampling  volume.  Most  spectral  analyses  have  been  based  on  the  Fast  Fourier  Transform 
(FFT).  However,  these  techniques  suffer  from  poor  spectral  resolution  and  a large 
variance  of  estimation.  Several  studies  have  shown  that  the  parametric  spectral  analysis 
can  provide  more  stable  spectral  estimation.  In  order  to  take  advantage  of  the  parametric 
spectral  analysis,  more  knowledge  is  required  on  the  Doppler  signal.  In  this  research,  the 
Doppler  signal  was  represented  as  a frequency-shifted,  moving  average  process.  A Data 
acquisition  system  was  designed  to  record  the  Doppler  signal.  Two  experiments  were 
conducted  for  the  Doppler  signal  from  the  middle  cerebral  artery  of  a normal  subject  and 
from  the  femoral  artery  of  a pig.  The  first  experiment  considered  the  design  problem 
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when  we  want  to  apply  ARSE  methods  to  the  Doppler  signal.  In  the  second  experiment, 
we  proposed  a blood  velocity  envelope  estimator  and  evaluated  the  performance. 

In  the  first  experiment,  five  analysis  tools  were  employed  to  solve  the  design 
problems.  Autocorrelation  and  the  Burg  methods  were  tested.  Three  sizes  of  analysis 
interval  (N)  were  test:  32,  64  and  128  points.  These  correspond  to  2.9,  5.8,  and  11.6 
msec.  The  autocorrelation  method  exhibited  more  stability  than  the  Burg  method.  The 
optimal  design  factors  were  found  using  the  autocorrelation  method  with  N=64  and  P=8. 

In  the  second  experiment,  a blood  velocity  envelope  estimation  technique  was 
proposed.  This  technique  consists  of  the  ARSE,  pole-tracing  algorithm,  and  a 
combination  of  median  and  three-point  Hamming  filter.  The  performance  of  this  technique 
under  various  conditions  was  evaluated  by  comparing  it  with  the  blood  flow  signal 
measured  by  an  electromagnetic  flowmeter.  The  autocorrelation  method  with  N=64  and 
P=8  yields  the  best  performance.  This  conforms  to  the  experimental  result  in  the  first 
experiment. 
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CHAPTER  1 


INTRODUCTION 
1.1.  Problem  Statement 

In  order  to  extract  quantitative  information  from  the  Doppler  signal,  many  studies 
have  been  conducted  on  the  theoretical  modeling  of  the  signal  [Brody  and  Meindl,  1974; 
Angelsen,  1980;  Garbini  et  al.,  1982;  Mo  and  Cobbold,  1986b],  In  these  models,  the 
characteristics  of  the  Doppler  signal  were  first  expressed  in  terms  of  an  autocorrelation 
function  (ACF).  The  Fourier  transform  of  the  ACF,  or  power  density  spectrum  (PSD), 
provides  quantitative  information  on  the  blood  velocity  distribution.  In  practice,  spectral 
analysis  has  been  accomplished  using  the  Fast  Fourier  Transform  (FFT)  algorithm,  which 
is  computationally  efficient.  However,  the  FFT  suffers  from  a large  variance  of  the 
estimate  due  to  the  implicit  windowing  of  the  data  and  poor  spectral  resolution.  These 
two  limitations  of  the  FFT  impose  serious  problems  when  analyzing  short  data  records. 
There  are  two  situations  that  give  rise  to  short  data  records:  when  available  data  are  brief 
in  duration;  and  when  the  data  are  time-varying  slowly  so  that  the  spectrum  may  be 
considered  constant  only  for  a short  analysis  interval.  The  spectral  analysis  of  the  Doppler 
signal  belongs  to  the  latter  case.  Namely,  the  Doppler  signal  is  a slowly  time-varying 
signal  so  that  the  Doppler  spectra  can  be  considered  constant  only  for  a short  analysis 
interval.  Due  to  the  drawbacks  of  TCD  signal  processing  and  the  complex  behavior  of  the 
blood  flow  in  human  arteries,  applications  of  the  TCD  have  been  limited  mainly  to 
qualitative  analyses.  Quantitative  analyses  have  employed  simple  heuristic  parameters. 

Alternative  spectral  estimation  methods  have  been  proposed  for  more  stable  spectral 
estimation  of  shorter  data  records  and  more  enhanced  spectral  resolution.  One  method  is 
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the  parametric  spectral  estimation.  The  parametric  spectral  estimation  method  estimates 
the  PSD  of  a stochastic  process  based  on  a linear  rational  model.  Therefore,  in  order  to 
take  advantage  of  the  parametric  spectral  analysis  methods,  more  knowledge  is  required 
from  the  Doppler  signal.  This  a priori  knowledge  of  the  model  is  the  key  to  successful 
spectral  estimation.  In  other  words,  the  structure  of  the  model  should  accurately  represent 
the  signal  generating  process. 

Among  the  many  parametric  spectral  estimation  methods,  the  autoregressive  (AR) 

method  has  attracted  great  interest  due  to  its  linear  solution  space  and  the  availability  of 

computationally  efficient  algorithms.  Kaluzynski  [1987]  reported  that  the  autoregressive 

spectra  were  more  stable  than  Fourier  spectra.  Kaluzynski  [1989]  and  Vaitkus  and 

Cobbold  [1988a,  1988b]  compared  Doppler  spectra  by  the  Fourier  method  and  the 

parametric  spectral  estimation  methods.  However,  previous  studies  were  purely 

experimental  and  incomplete.  Most  of  the  parametric  spectral  estimation  studies  were 

performed  on  the  direct  Doppler  signal  without  any  reference  to  rigorous  mathematical 

signal  models.  For  example,  the  study  of  Kaluzynski  [1987]  was  purely  experimental  and 

provided  no  theoretical  rationale  for  the  AR  methods.  Also,  he  did  not  address  the 

question  of  accuracy  or  the  effect  of  model  order  as  well  as  the  size  of  the  analysis 

interval.  His  1989  paper  [Kaluzynski,  1989]  referred  to  the  Wold  theorem,  which  allows 

the  application  of  an  AR  model  of  sufficiently  high  order  to  any  signal  that  can  be 

presented  by  a linear  model.  However,  he  did  not  provide  any  evidence  supporting  the 

linear  model  that  generates  the  Doppler  signal.  In  their  comprehensive  investigation  of 

parametric  spectral  analysis  of  the  Doppler  signal,  Vaitkus  and  Cobbold  [1988a]  intuitively 

suggested  that  the  Doppler  signal  could  be  considered  to  arise,  in  part,  from  a moving 

average  process.  Although  their  rationale  is  intuitive,  it  is  worth  quoting 

The  Doppler  signal  arises  from  the  scattered  contribution  from  all  cells  within  the 
beam.  As  time  progresses,  the  movement  of  new  cells  into  the  beam  volume  and 
old  cells  out  causes  the  distribution  of  cells  to  change,  thereby  resulting  in  the 
Doppler  signal  becoming  decorrelated  until,  after  a time  equal  to  the  transit  time 
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of  the  cells  through  the  beam,  there  is  an  entirely  new  distribution  of  scatterers. 

After  this,  the  signal  is  completely  uncorrelated  with  that  occurring  earlier 
[Vaitkus  and  Cobbold,  1988a,  pp.  664-665], 

Yet,  this  rationale  was  not  sufficiently  developed  into  a rigorous  mathematical  model.  In 
order  to  take  advantage  of  the  parametric  spectral  estimation  method,  we  must  start  with 
an  accurate  structure  of  a model  that  represents  a signal  generating  process.  The  lack  of 
mathematical  rigor  in  describing  the  modeling  of  the  Doppler  signal  motivated  this 
research  effort. 


1.2.  Background 

An  Ultrasound  Doppler  system  (UDS)  has  been  used  to  measure  blood  velocity  inside 
the  human  body.  A Transcranial  Doppler  (TCD)  device  specializes  in  measuring  the  blood 
velocity  inside  major  arteries  of  the  human  brain,  or  cerebral  arteries.  Aaslid  et  al.  [1982] 
first  reported  the  applications  of  the  TCD  in  the  human  brain,  using  ultrasound  with  a 
frequency  of  1-2  MHz  (versus  the  5-10  MHz  of  other  clinically  used  ultrasound  Dopplers) 
to  decrease  the  attenuation  caused  by  bone  and  soft  tissue.  Since  the  TCD  is  noninvasive 
and  real-time,  it  has  advantages  over  other  blood  flow  measurement  modalities  such  as  an 
electromagnetic  flowmeter.  Accurate  knowledge  of  the  pulsatile  blood  flow  velocity 
provides  invaluable  information  for  diagnosis  of  cerebrovascular  disease  [Aaslid  et  al., 
1984]  and  for  monitoring  of  blood  supply  in  human  brain  during  anesthesia  [Otis  and 
Ringelstein,  1990;  Halsey,  1989], 

Figure  1-1  is  a functional  block  diagram  of  the  TCD.  The  UDS  operates  in  two  types 
of  operating  modes:  pulsed  wave  (PW)  and  continuous  wave  (CW)  modes.  The  TCD 
operates  only  in  the  PW  mode.  The  PW  mode  enables  the  measurement  of  the  blood 
velocity  in  a specific  segment  of  an  artery,  which  is  called  the  sample  volume.  This 
capability  is  called  target  ranging  in  the  field  of  RADAR.  Target  ranging  is  not  possible  in 
the  CW  mode  because  the  received  signal  may  contain  the  ultrasound  scattered  from  the 
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Fig.  1-1  A block  diagram  of  Transcranial  Doppler 
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blood  in  various  arteries  that  lies  within  the  path  of  propagation.  The  main  reasons  why 
the  TCD  requires  target  ranging  are 

1)  There  are  numerous  arteries  inside  the  human  brain. 

2)  The  ultrasound  wave  can  penetrate  human  skull  only  through  limited  regions. 
Therefore,  the  capability  of  target  ranging  is  essential  to  pinpoint  a specific  artery  among 
numerous  cerebral  arteries.  Aaslid  et  al.  [1982]  showed  that  most  major  cerebral  arteries 
can  be  accessed  through  the  temporal  region,  since  the  bone  of  the  temporal  region  is 
thinner  than  others. 

In  PW  mode,  the  transducer  transmits  a pulse-modulated  ultrasound  wave  into  the 
Ml  section  of  the  MCA.  Subsequently,  blood  in  the  sample  volume  backscatters  the 
ultrasound  wave.  It  is  well  known  that  the  red  blood  cells  (RBCs)  are  the  main  scatterers 
of  the  ultrasound  wave  [Reid  et  al.,  1969],  If  the  red  blood  cells  are  moving,  the 
backscattered  ultrasound  is  frequency-shifted  by  the  Doppler  effect.  The  frequency  shift, 
fd,  by  a single  scatterer  can  be  written  as 


2vf0  cos<J) 
c 


(1-1) 


where 

v:  scatterer  velocity 
c:  speed  of  ultrasound 
f0:  transmitted  frequency 

<j):  angle  between  transmitted  beam  and  the  direction  of  the  moving  scatterer. 
Subsequently,  the  backscattered  ultrasound  is  received  by  the  same  transducer.  By 
controlling  the  time  between  the  transmission  and  the  reception,  distance  from  the 
transducer  and  the  sample  volume  can  be  varied.  The  received  signal  is  passed  through  a 
frequency  demodulator,  thereby  leaving  only  the  Doppler  frequency  shift  components. 
The  output  of  the  demodulator  is  called  the  Doppler  signal. 
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The  PW  Doppler  signal  from  a human  artery  has  been  modeled  as  a nonstationary, 
discrete  stochastic  process.  Since  the  position  and  the  reflection  properties  of  a red  blood 
cell  can  only  be  characterized  in  probabilistic  terms,  the  Doppler  signal  has  been 
represented  as  a stochastic  process.  Because  of  the  pulsatile  nature  of  the  blood  flow,  the 
Doppler  signal  is  nonstationary.  However,  the  Doppler  signal  has  been  considered  to  be  a 
wide-sense  stationary  process  so  that  it  can  be  assumed  to  be  stationary  within  a small 
analysis  interval  of  10  msec  [Mo  and  Cobbold,  1986a], 

The  power  spectral  density  (PSD)  of  the  Doppler  signal  is  mainly  affected  by  blood 
velocity  distribution  within  a sample  volume.  Therefore,  the  bandwidth  of  the  Doppler 
signal  may  be  broad  or  narrow  depending  on  the  velocity  distribution.  The  PSD  of  the 
Doppler  signal  is  also  affected  by  other  factors  that  results  in  broadening  the  bandwidth 
[Newhouse  et  al.,  1976;  Morin  et  al.,  1988],  These  factors,  considered  unwanted  noise, 
are  called  the  spectral  broadening  factors.  One  of  the  spectral  broadening  factors,  called 
the  transit  time  effect,  is  due  to  the  finite  duration  of  the  moving  red  blood  cells  in  a 
sample  volume.  Geometric  spectral  broadening  is  caused  by  the  dimension  of  the 
transducers.  These  spectral  broadening  factors  are  inherent  phenomena  of  the  ultrasound 
Doppler  instruments. 

In  order  to  assess  the  velocity  information,  we  need  to  process  the  Doppler  signal. 
There  are  three  distinctive  processing  methods  of  that  are  usually  applied  onto  the 
Doppler  signal 

1)  audio  analysis, 

2)  time-domain  analysis, 

3)  and  spectral  analysis. 

The  most  primitive  of  the  methods  is  to  listen  to  the  Doppler  signal,  since  the  frequency  of 
the  Doppler  signal  is  within  the  audible  range  of  human  ears  [Beach  and  Nix,  1987], 
Although  the  audio  analysis  is  simple  and  economical,  it  is  largely  qualitative  and  requires 
expert  knowledge  for  interpretation.  Some  time-domain  analysis  techniques  are  based  on 
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the  definition  of  the  instantaneous  frequency  or  mean  frequency  for  the  stationary  signal. 
The  instantaneous  frequency  is  defined  as  the  mean  of  time-frequency  distribution  at  a 
specific  time.  The  mean  frequency  is  defined  as  the  mean  of  frequency  distribution  during 
an  interval.  The  physical  meanings  of  both  instantaneous  and  mean  frequency  is  the  mean 
of  velocity  distribution  within  a sample  volume. 

The  short-time  spectral  analysis  resorts  to  the  explicit  reference  to  the  Fourier 
transform  of  the  Doppler  signal  (or  the  autocorrelation  of  the  Doppler  signal).  This 
results  in  a power  spectral  density  (PSD)  of  the  Doppler  signal  during  an  analysis  interval. 
The  physical  meaning  of  the  Doppler  PSD  is  the  velocity  distribution  with  in  a sample 
volume. 

From  the  viewpoint  of  medical  applications,  spectral  analysis  falls  into  two  distinctive 
research  areas.  The  first  area  concentrates  on  estimating  the  spectrum  of  the  Doppler 
signal  itself  [Vaitkus  and  Cobbold,  1988a,  1988b;  Kaluzynski,  1987,  1989].  Another  issue 
of  this  area  is  to  develop  simulation  models  in  order  to  verify  the  proposed  signal  analysis 
methods  [Mo  and  Cobbold,  1986a;  Jones  and  Giddens,  1990],  Many  spectral  analysis 
techniques  such  as  the  periodogram,  the  Blackman-Tukey  method,  parametric  spectral 
analysis,  and  the  Wigner-Ville  distribution,  are  assessed  for  the  velocity  estimation.  Figure 
1-2  is  a typical  gray-level  display  of  the  PSD  of  the  Doppler  signal.  The  horizontal  and 
the  vertical  axis  indicate  the  time  and  frequency  of  the  signal,  respectively.  The  PSD, 
found  by  the  Fast  Fourier  Transform  (FFT),  is  mapped  to  a gray-level  of  each  pixel  on  the 
display.  The  main  application  of  this  area  is  the  diagnosis  of  arterial  deformation,  because 
the  PSD  of  the  Doppler  signal  found  around  the  peak  systole  have  been  shown  to  give 
quantitative  information  for  the  assessment  of  carotid  arterial  disease  [Kalman  et  al., 

1985], 

The  other  area  deals  with  the  derivation  of  the  flow  velocity  waveform  [D'Alessio, 
1985;  Kitney  et  al.,  1986;  Mo  et  al.,  1988],  The  spectral  analysis  of  the  Doppler  signal 
yields  the  velocity  of  flow  (cm/sec)  rather  than  the  volume  flow  (cmVsec).  However,  most 
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Fig.  1-2  Typical  gray-level  display  of  the  TCD 
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clinicians  are  interested  in  the  volume  flow  rather  than  the  velocity.  For  example, 
monitoring  the  blood  flow  to  the  human  brain  may  be  of  great  importance  during  general 
anesthesia.  The  mean  of  the  Doppler  spectrum  may  be  the  first  candidate  for  the 
derivation  of  the  blood  flow.  However,  other  parameters  such  as  mode  and  maximum 
frequency  can  be  used.  The  most  popular  method  in  this  area  is  the  maximum  velocity 
envelope  (MVE)  estimation  [Thompson  et  al.,  1986;  Marlcwalder  et  al.,  1984],  The  MVE 
can  be  thought  of  as  the  velocity  of  the  fastest  moving  blood  cells  within  the  sample 
volume.  There  are  several  ways  to  derive  the  MVE,  but  most  of  them  are  based  on  the 
FFT  [D'Alessio,  1985;  Mo  et  al.,  1988], 

1.3.  Research  Goals  and  Approach 

The  goals  of  this  research  are 

1)  to  develop  a theoretical  parametric  model  of  the  PW  Doppler  signal, 

2)  to  determine  the  optimal  parameters  for  autoregressive  spectral  estimation 
(ARSE),  and 

3)  to  develop  a velocity  envelope  estimation  scheme  based  on  the  ARSE  methods. 

If  one  is  successful  in  developing  a parametric  model  for  the  behavior  of  some  signals,  the 
model  can  enhance  the  abilities  for  various  applications  such  as  prediction,  data 
compression,  and  spectral  estimation  with  high  resolution.  In  speech,  for  example,  it  is 
known  that  the  speech  generating  mechanism  can  be  modeled  as  an  all-pole  linear  system 
[Markel  and  Gray,  1976],  This  all-pole  model  has  been  used  successfully  for  applications 
such  as  speech  synthesis;  speech  recognition;  and  speech  data  compression.  Therefore, 
for  our  Doppler  signal  processing,  a parametric  representation  of  the  signal  is  the  essential 
first  step  for  the  applications  of  parametric  spectral  analysis. 

The  parametric  model  of  the  Doppler  signal  is  based  on  the  model  proposed  by  Jones 
and  Giddens  [1990],  They  developed  a model  to  simulate  the  PW  Doppler  signal  in  order 
to  assess  the  performance  of  spectral  estimation  such  as  the  Fourier  method  and 
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autoregressive  method.  In  this  study,  we  project  the  model  from  a different  point  of  view. 
We  can  use  the  simulation  model  as  the  basic  tool  to  design  a parametric  spectral 
estimation.  That  is,  the  simulation  model  will  provide  physical  meaning  of  the  model 
order  and  the  system  parameters  of  the  spectral  estimator. 

In  order  to  determine  optimal  parameters  for  the  ARSE  model,  the  design  factors  for 
the  autoregressive  spectral  estimation  are  determined  using  five  analysis  procedures:  1) 
final  prediction  error  (FPE),  2)  criterion  autoregression  transfer  (CAT),  3)  Akaike 
information  criterion  (AIC),  4)  spectral  distance  measure  (SDM),  and  5)  spectral  flatness 
measure  (SFM).  As  the  first  step,  the  rationale  for  the  autoregressive  modeling  of  the 
Doppler  signal  is  discussed.  Two  algorithms  (autocorrelation  and  Burg  methods)  are 
tested.  Three  sizes  of  the  analysis  interval  (N)  were  tested:  32,  64  and  128  points.  These 
correspond  to  2.9,  5.8,  and  11.6  [msec]. 

Finally,  a blood  velocity  envelope  estimation  technique  is  proposed.  This  technique 
consists  of  the  ARSE,  pole-tracing  algorithm,  and  a combination  of  a median  and  three- 
point  Hamming  filter.  The  performance  of  this  technique  under  various  conditions  was 
evaluated  by  comparing  it  with  the  blood  flow  signal  measured  by  an  electromagnetic 
flowmeter.  Time-domain  analysis  technique  (the  complex  autocorrelation  method)  was 
also  included  for  evaluation.  The  evaluation  involved  a linear  regression  of  two 
waveforms.  The  variance  of  the  error  was  used  for  the  quantitative  evaluation. 

1.4.  Contributions 

The  proposed  research  has  the  following  potential  benefits 

1)  The  development  of  a parametric  spectral  analysis  of  the  Doppler  signal. 

2)  A methodology  for  the  optimal  design  of  the  TCD. 

3)  A more  rigorous  mathematical  description  of  the  underlying  process  generating 


the  Doppler  signal 
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1.5.  Description  of  Chapters 

Chapter  2 presents  a review  of  the  research  areas  associated  with  ultrasound  Doppler 
instruments  in  the  medical  domain.  The  research  is  divided  into  three  areas:  Doppler 
signal  modeling,  Doppler  signal  processing,  and  the  medical  applications  of  the  ultrasound 
Doppler  instruments.  Considerable  attention  was  given  to  the  time-domain  analysis  and 
maximum  velocity  envelope  estimation  of  the  Doppler  signal. 

Chapter  3 presents  the  general  theory  behind  parametric  spectral  estimation, 
especially  the  autoregressive  spectral  estimation  method  (ARSE).  The  linear  prediction 
formation  of  the  ARSE  method  will  be  presented  in  detail.  Other  formulations  which 
result  in  the  ARSE  will  also  be  reviewed.  Implementation  algorithms  to  estimate  the  AR 
model  parameters  with  limited  data  records  are  presented  next. 

In  Chapter  4,  a parametric  model  of  the  Doppler  signal  is  introduced.  The 
assumptions  required  to  perform  the  modeling  are  discussed.  Subsequently,  a 
mathematical  representation  is  derived  in  terms  of  the  characteristics  of  both  the  UDS 
parameters  and  the  blood  flow. 

In  Chapter  5,  the  hardware  and  software  implementation  for  data  collection  and 
recording  of  the  Doppler  signal  is  described.  The  implemented  data  structures  are 
discussed.  The  control  flow  of  the  software  is  also  discussed. 

In  Chapter  6,  the  design  factors  of  the  ARSE  are  considered.  First,  rationale  for  the 
ARSE  of  the  Doppler  signal  is  provided.  Several  analysis  procedures  are  employed  in 
order  to  determine  the  optimal  values  of  the  model  order  and  the  analysis  interval. 

A velocity  envelope  estimator  based  on  the  ARSE  method  is  proposed  in  Chapter  7. 
The  performance  of  our  technique  is  compared  with  the  blood  flow  signal  measure  by  an 
electromagnetic  flowmeter. 

Chapter  8 summarizes  our  conclusions  in  this  study  and  gives  suggestions  for  future 


research. 


CHAPTER  2 


REVIEW  OF  ULTRASOUND  DOPPLER  TECHNOLOGY 
IN  THE  MEDICAL  DOMAIN 

2.1.  Introduction 

In  this  chapter,  research  areas  associated  with  ultrasound  Doppler  system  (UDS)  in 
the  medical  domain  are  reviewed.  Early  research  is  reviewed  which  attempted  to  apply 
ultrasound  Doppler  technology  in  humans.  Our  discussion  is  focused  in  studies  which  fall 
into  three  groups  as  shown  in  Fig.  2-1,  namely,  Doppler  signal  modeling,  Doppler  signal 
processing,  and  medical  applications.  Doppler  signal  modeling  attempts  to  describe  the 
Doppler  signal  in  terms  of  the  physical  characteristics  of  the  UDS  and  the  hemodynamic 
properties  of  the  blood  in  a sample  volume.  This  includes  the  mathematical  modeling  of 
each  of  the  blocks  of  Fig.  1-1  that  generate  the  Doppler  signal,  mainly  the  transmitter, 
transmitting  transducer,  scattering  process,  receiving  transducer,  and  the  demodulator. 

The  Doppler  signal  processing  attempts  to  obtain  a more  convenient  representation  of  the 
information  carried  by  the  Doppler  signal.  Parameters  measured  include  the  mean  blood 
flow,  the  velocity  profile  and  the  existence  of  turbulent  blood  flow.  The  various  clinical 
applications  of  the  UDS  impose  different  design  criteria  to  both  the  physical  characteristics 
of  the  UDS  and  the  Doppler  signal  processing.  For  example,  the  information  in  the 
velocity  profile  requires  a small  sampling  volume  and  is  vital  to  the  diagnosis  of  a stenosis 
(partial  and  complete  blockage  in  an  artery).  In  order  to  decrease  the  size  of  the  sample 
volume,  the  UDS  should  be  operated  in  the  pulse-wave  mode  with  a short  pulse  width. 
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Fig.  2-1  Three  research  areas  in  Ultrasound  Doppler  Technology 
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2 2 Early  Research  in  Ultrasound  Doppler  Techniques 
The  application  of  the  ultrasound  Doppler  technology  to  humans  was  first  originated 
by  Satomura  [Satomura,  1959;  Kaneko,  1986],  Before  Satomura,  ultrasound  had  been 
used  for  30  years  in  the  military  for  the  detection  of  submarines.  The  initial  purpose  of 
Satomura's  study  was  to  assess  the  movement  of  a human  organ  such  as  the  heart  motion 
and  the  pulsations  of  blood  vessel  walls.  A continuous-wave  (CW)  ultrasound  with  a 
known  frequency  was  insonated  in  a moving  organ.  By  measuring  the  Doppler  frequency 
shift  of  the  returned  ultrasound  wave,  Satomura  could  determine  the  movement.  During 
the  measurement  of  the  vibration  of  the  blood  vessel  wall,  he  also  found  high  frequency 
Doppler  noise.  He  confirmed  that  the  noise  was  not  produced  by  the  vessel  walls  but  by 
the  blood  flow.  This  led  to  the  development  of  the  ultrasound  Doppler  system  for 
measuring  blood  velocity.  He  speculated  that  the  source  of  the  noise  was  the  turbulent 
blood  flow  in  the  blood  stream.  Later,  Kato  et  al.  [1962]  established  that  the  Doppler 
noise  was  caused  by  the  blood  cells  and  not  by  the  turbulent  blood  flow.  In  the  United 
States,  Flanklin  et  al.  [1961]  first  introduced  a paper  describing  a CW  ultrasound  blood 
flowmeter.  They  used  to  detect  blood  flow  in  animals. 

Clinicians  began  also  to  ask  more  sophisticated  diagnostic  questions.  As  their 
understanding  of  the  dynamic  blood  flow  expanded,  they  desired  clinically  important 
measurements  of  the  site  of  the  occlusions,  volume-flow  rates,  velocity  profiles,  and 
volume  profiles.  One  of  the  principal  limitations  of  the  early  CW  Dopplers  was  their 
inability  to  detect  the  direction  of  the  flow  with  respect  to  the  transducer.  In  1966, 
McCleod  demonstrated  the  first  direction-sensing  CW  Doppler  device.  He  used  a version 
of  the  quadrature  phase  detection  common  to  single  sideband  communications  to  separate 
the  upper  and  lower  Doppler  sidebands  from  the  carrier.  Another  advancement  was  the 
pulsed-wave  (PW)  ultrasound  Doppler.  The  PW  Doppler  had  an  important  advantage 
over  the  CW  ultrasound  Doppler  in  target  ranging.  Baker  [1970]  worked  on  the 
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implementation  of  a PW  ultrasound  Doppler.  Using  the  device,  a flow  velocity  profile  was 
mapped  over  the  period  of  pulsatile  flow. 

2.3.  Doppler  Signal  Modeling 

Doppler  signal  modeling  expounds  on  the  mathematical  modeling  of  each  block  of 
Fig.  1-1.  Significant  theoretical  research  in  the  ultrasound  Doppler  focused  on  the 
scattering  process.  The  scattering  process  describes  the  interaction  between  the  blood  and 
the  ultrasound  wave.  As  new  experimental  findings  on  the  blood  scattering  were 
discovered  [Reid  et  al.,  1969;  Shung  et  al.,  1976],  new  models  of  the  scattering  process 
were  developed  [Brody  and  Meindl,  1974;  Angelson,  1980;  Mo  and  Cobbold,  1986b].  In 
order  to  understand  the  interaction  of  ultrasound  and  blood,  researchers  performed  many 
experiments  varying  the  transmitting  carrier  frequency,  hematocrit,  laminar  and  turbulent 
blood  flows. 

The  first  important  experiment  on  the  scattering  of  the  ultrasound  by  the  human  blood 
was  reported  by  Reid  et  al.  [1969],  Their  findings  were  compatible  with  the  assumption 
that  the  scattering  of  the  ultrasound  by  blood  obeys  the  Rayleigh  scattering  law.  That  is,  if 
the  particle  size  is  much  less  than  the  wavelength  of  the  ultrasound  and  if  the  scattering  is 
first  order  and  isotropic,  then  the  scattering  energy  will  increase  as  the  fourth  power  of  the 
incident  wave  frequency.  These  experimental  results  played  an  essential  role  providing 
reasonable  assumptions  for  the  theoretical  modeling  of  the  scattering  process. 

The  paper  by  Brody  and  Meindl  [1974]  was  the  first  general  and  comprehensive 
theoretical  analysis  of  a continuous-wave  ultrasound  Doppler  signal  for  steady,  laminar 
blood  flow.  They  identified  basic  assumptions  required  to  develop  a general  model.  Their 
formulation  was  greatly  influenced  by  the  statistical  studies  on  reverberation  in  Sonar 
[Faure,  1964;  Van  Trees,  1965;  Middleton,  1967],  The  reverberation  is  considered  as 
noise  produced  by  scatterers  in  a random  medium.  The  experimental  results  by  Reid  et  al. 
[1969]  provided  reasonable  ground  to  develop  a model  of  blood  as  a random  medium  in  a 
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way  similar  to  the  reverberation  model.  Brody  and  Meindl  first  assumed  that  each  red 
blood  cell  (RBC)  is  an  identical  independent  point  scatterer.  In  addition,  the  amplitude 
and  phase  shift  of  the  scattered  radiation  by  an  RBC  were  assumed  to  be  random  [Faure, 
1964;  Van  Trees,  1965],  Using  these  and  other  assumptions,  they  developed  a stochastic 
model  which  characterized  the  back-scattered  ultrasound  as  a Gaussian  random  process. 
An  expression  for  the  autocovariance  function  of  the  Doppler  signal  was  also  obtained. 

In  their  study  on  the  scattering  of  the  ultrasound  by  blood,  Shung  et  al.  [1976]  found 
that  the  scattering  coefficient  increased  along  with  the  hematocrit  (percentage  of  the  blood 
cells  in  a unit  volume)  until  it  reached  a maximum  around  26  percent  and  then  decreased 
as  the  hematocrit  increased.  The  normal  hematocrit  of  a human  is  45  percent,  and  the 
average  distance  between  two  RBCs  is  about  10  percent  of  its  diameter.  Thus,  the  RBCs 
can  not  be  considered  independent  of  each  other.  They  speculated  that  beyond  a certain 
point,  the  scattering  process  became  more  complicated  so  that  a single  scatterer  would  not 
suffice  to  describe  the  scattering  mechanism.  The  complication  was  speculated  to  be 
caused  by  either  multiple  scattering  or  the  disappearance  of  the  random  nature  of  the 
motion  of  scatterers. 

The  experimental  results  of  Shung  et  al.  [1976]  led  to  a serious  problem  in  applying 
the  assumption  of  independence  in  the  model  by  Brody  and  Meindl.  Angelson  [1980] 
developed  a model  with  the  assumptions  that  the  scattering  was  first-order  in  nature  but 
that  each  scatterer  was  dependent  in  position.  However,  he  considered  normal  blood  as  a 
continuum  composed  of  base  cells  so  that  beyond  the  boundary  of  a cell,  RBCs  were  no 
longer  dependent.  The  main  theme  of  Angleson's  formulation  was  that  we  could  still 
consider  these  independent  base  cells  as  point  scatterers  compared  to  the  wavelength  of 
the  insonated  ultrasound.  In  his  model,  the  scattering  arises  from  fluctuations  in  the  mass 
density  and  compressibility  of  the  blood,  which  is  caused  by  the  fluctuation  in  RBC 
concentration. 
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The  paper  by  Baker  [1970]  addressed  the  instrumentation  of  the  PW  ultrasound 
Doppler.  Garbini  et  al.  [1982]  developed  a stochastic  model  of  the  PW  ultrasound 
Doppler  for  homogeneous  turbulence  in  steady  flow.  One  of  their  contributions  was  to 
mathematically  clarify  the  effects  of  the  transit  time  effect  which  they  called  the  Doppler 
ambiguity  process. 

In  1986,  Mo  and  Cobbold  published  a general  model  of  the  CW  ultrasound  Doppler. 
Angelson  [1980]  included  the  effects  of  an  interacting  (or  correlated)  scatterer.  Mo  and 
Cobbold  [1986]  went  one  step  further  by  including  the  effects  of  blood  aggregation.  That 
is,  they  assumed  that  the  scatters  were  not  identical. 

In  1990,  Jones  and  Giddens  proposed  a model  for  the  PW  Doppler  ultrasound.  This 
model  was  used  to  simulate  a Doppler  signal.  An  important  component  of  the  model  was 
a sliding  window  used  to  correlate  normally  distributed  parameters.  This  window  proved 
to  be  related  to  the  shape  of  a sample  volume  and  transmitted  pulse. 

In  summary,  the  details  of  the  interaction  between  the  ultrasound  and  the  human 
physiology  are  not  clear.  The  main  problems  is  the  difficulty  describing  the  scattering 
process.  The  information  which  the  ultrasound  Doppler  signal  conveys  is  still  a mystery  of 
sorts.  Furthermore,  complex  hemodynamic  behavior  hinders  evaluation  of  the 
measurements. 


2,4.  Doppler  Signal  Processing 

The  main  purpose  of  Doppler  signal  processing  is  to  obtain  blood  velocity 
information.  The  representation  may  be  qualitative  or  quantitative.  As  in  Fig.  2-2,  there 
are  three  popular  ways  to  represent  velocity  information  involving  audio  analysis,  time- 
domain  analysis,  and  spectral  analysis.  Audio  analysis  is  the  oldest  method  requiring 
spectral  evaluation  by  human  scorers.  This  analysis  is  mainly  qualitative.  Time-domain 
analysis  attempts  to  calculate  mean  and  instantaneous  blood  velocity  as  a function  of  time. 
Numerous  laboratory  studies  have  shown  that  the  time-domain  analysis  depends  on  the 
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Fig.  2-2  Doppler  signal  processing  methods 
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velocity  profile  which  exhibited  variation  among  different  arteries  and  among  different 
patients.  The  Fast  Fourier  Transform  (FFT)  algorithm  [Cooley  and  Tukey,  1965]  made 
spectral  analysis  of  the  Doppler  signal  possible.  The  following  sections  describe  each 
method. 

2.4.1.  Audio  Analysis 

Audio  analysis  relies  on  spectral  discrimination  by  human  listeners.  It  has  been  used 
as  a qualitative  tool  to  identify  a blood  vessel  and  the  presence  of  stenosis  [Beach  and  Nix, 
1987],  For  example,  a trained  examiner  is  able  to  identify  the  three  branches  of  the  carotid 
artery  found  in  the  neck  by  their  audible  characteristics.  The  internal  carotid  artery, 
supplying  a low  resistance  vascular  bed,  has  almost  continuous  forward  flow  throughout  a 
heart  beat.  The  external  carotid  artery,  supplying  a high  resistance  vascular  bed,  has 
multiphase  qualities.  The  common  carotid  artery,  which  supplies  both  the  internal  and 
external  arteries,  has  audio  characteristics  that  fall  between  the  internal  and  external 
arteries.  The  examiner  can  also  assess  the  artery  and  its  branches  for  significant  increases 
in  frequency,  or  pitch,  of  the  Doppler  signal.  In  addition  to  the  audible  high  pitch 
associated  with  high  blood  velocity,  the  examiner  can  differentiate  the  "hiss"  associated 
with  spectral  broadening  from  the  rapidly  fluctuating  "gurgling"  sound  associated  with 
flow  disturbance.  The  frequency  variations  followed  by  a disturbed  signal  are  associated 
with  a stenotic  area.  Nix  et  al.  [1979]  reported  that  the  examiner  can  identify  at  least  70% 
of  those  patients  with  a high-grade  stenosis  condition  equal  to  50%  diameter  reduction  or 
greater. 

Audio  analysis  is  the  most  economic  way  to  assess  the  Doppler  signal.  However,  it 


demands  expert  listeners  and  scorers. 


20 


2.4.2.  Time-Domain  Analysis 

In  the  early  research  days  of  Doppler  signal  processing,  inexpensive  real-time  spectral 
Fourier  analysis  instruments  were  not  available.  Signal  processing  relied  on  time-domain 
analysis  to  obtain  the  mean  blood  velocity  and  the  variance  or  bandwidth  of  the  Doppler 
spectrum  which  are  clinically  significant.  Research  papers  on  the  time-domain  Doppler 
signal  processing  emphasized  the  simplicity  of  the  implementation,  the  possibility  of  real- 
time signal  processing,  the  performance  under  various  clinical  conditions,  and  the  validity 
of  the  blood  velocity  estimation  in  various  arteries.  We  will  first  review  the  basic 
principles  used  in  digital  time-domain  analysis  of  the  Doppler  signal.  Various 
implementation  techniques  are  discussed,  followed  by  a performance  assessment  of  each 
implementation  option. 

2.4,2. 1.  Principle 

The  received  signal  in  a CW  UDS  can  be  expressed  mathematically  as 

r(t)  = Re[z(t)  exp(jco  0t)]  (2-1) 

where  co0  is  the  transmitted  angular  frequency,  and  z(t)  is  the  complex  envelope.  The 
complex  envelope  contains  both  amplitude  and  phase  information  of  the  received  signal. 
For  stationary  blood  flow,  this  envelope  can  be  modeled  as  a stationary  complex  Gaussian 
process,  the  rationale  of  which  is  based  on  the  assumption  that  the  signal  is  composed  of 
the  sum  of  uncorrelated  contributions  from  a large  number  of  scatterers  [Brody  and 
Meindl,  1974;  Angelsen,  1980],  The  complex  envelope  can  be  expressed  mathematically 
as 

z(t)  = x(t)  + jy(t)  (2-2) 

where  x(t)  and  y(t)  are  called  the  quadrature  components  of  the  Doppler  signal.  The 
complex  Gaussian  process  in  Eq.  (2-2)  can  be  expressed  in  terms  of  amplitude  a(t)  and 
phase  0(t)  such  that 
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z(t)  = x(t)  + jy(t)  = a(t)  exp(j0(t)) 


(2-3) 


where 


a(t)  = Vx2( t)  + y2(t),  and  0(t)  = tan1  ^ 


x(t) 


(2-4) 


In  the  PW  mode,  the  signal  is  a sampled  version  of  the  continuous  signal  with  the  sampling 


frequency  equal  to  the  pulse  repetition  frequency  (PRF).  In  this  case,  we  can  regenerate 
the  analog  signal  with  a low-pass  filter,  or  we  can  implement  a digital  system  which 
processes  the  digital  Doppler  signal  directly. 

The  instantaneous  angular  frequency,  co^t),  of  the  process  is  defined  as  the  time 
derivative  of  the  phase  0(t),  that  is, 


(2-5) 


_ x(t)y*(t)-x,(t)y(t) 
x2  (t)  + y2  (t) 


where  x'(t)  and  y'(t)  are  the  time  derivatives  of  x(t)  and  y(t),  respectively. 
The  mean  angular  frequency,  com,  of  the  process  is  defined  as 
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where  S^co)  is  the  power  spectral  density  (PSD)  of  z(t),  and  1^(0)  is  the  autocorrelation 
function  (ACF)  of  x(t)  with  zero  lag.  The  PSD  is  defined  as  the  Fourier  transform  of  the 
ACF,  that  is 


22 


S„(o)  = | R^x)  exp(-j  ©x)dx  (2-7) 

-oo 

and  the  ACF  of  z(t)  is  defined  as 

Rzz(t)  =<  z*(t)z(t  + x)  > (2-8) 

where  < • > denotes  ensemble  averaging,  and  * indicates  complex  conjugate.  Using  Eq. 
(2-3),  the  ACF  of  z(t)  can  be  expressed  in  terms  of  the  ACFs  of  its  quadrature 
components  such  that 

R„(t)  = R„(t)  + R„(t)+j  [R„(t)-R„(t)].  (2-9) 

Then,  we  can  express  com  as 


„ = -R'y.<°)  (2.10' 

” j *»(«)  R„(0)+  R„(0) 

where  the  last  equality  is  obtained  from  the  facts  that  R^  and  R^  are  even  functions,  while 
R*y  and  R^  are  odd,  and  that  derivative  of  an  even  function  is  odd,  and  vice  versa.  The 
above  equation  can  expressed  in  terms  of  x(t)  and  y(t) 


< x(t)y'(t)  > - < x'(t)y(t)  > 
< x2(t)  + y2(t)  > 


(2-11) 


Note  that  com  is  obtained  by  averaging  the  numerator  and  denominator  in  the  expression 
for  C0j(t).  In  practice,  we  assume  that  the  Doppler  signal  is  ergodic,  and  the  ensemble 
average  is  replaced  by  the  time  average  and  we  can  write  Eq.  (2-11)  such  that 
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J(x(t)  y'(t)-x'(t)  y(t))dt 
-00 

J(x!(t)  + yI(t))dt 

-oo 


(2-12) 


For  the  complex  envelope  to  be  stationary,  the  quadrature  components  must  have  the 
following  properties  [Papoulis,  1981] 

r„«=r„«,  *,(')■-*,(')■  e-13) 

Under  these  conditions,  the  autocorrelation  function  of  z(t)  then  becomes 

R=W=2(R„W+jRwW).  (2-14) 

Using  this  expression,  the  mean  angular  frequency  is  expressed  as 
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_<x'(t)  y(t)  > 
<x2(t)>  ' 


(2-15) 


The  definition  of  the  instantaneous  frequency  is  only  valid  for  monocomponent  signals. 
When  the  signal  is  multicomponent  or  a broad  band  stochastic  signal,  the  instantaneous 
frequency  should  be  a good  estimation  of  the  mean  value  of  the  time-frequency 
distribution  at  a specific  time.  Angelsen  [1981]  showed  that  for  a deterministic  signal,  the 
instantaneous  frequency  is  equal  to  the  mean  value  of  the  frequency  distribution  plus  an 
oscillatory  term  that  contains  a time  derivative  of  the  amplitude  modulation.  Therefore,  if 
the  signal  is  narrowband  or  if  we  average  the  instantaneous  frequency  over  a short  period 
of  time,  the  oscillatory  term  can  be  neglected. 

The  variance  of  a spectrum,  a2 , is  defined  as 
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In  terms  of  the  ACF,  the  variance  can  be  written  as 
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'R'„(0)Yr"„(0) 
R„(0)  J R„(0)  ' 


(2-17) 


2.4. 2.2,  Digital  implementation 

There  are  various  methods  to  implement  a mean  frequency  estimator  using  digital 
techniques.  One  implementation  group  starts  with  a definition  of  the  mean  frequency  in 
terms  of  the  quadrature  components  as  in  Eq.  (2-12).  Two  methods  are  described:  the 
phase  detector  and  the  instantaneous  frequency  estimator.  Another  implementation  is 
based  on  the  expression  of  the  mean  frequency  in  terms  of  the  ACF.  This  implementation 
is  called  the  complex  autocorrelator,  or  simply  the  autocorrelator.  Finally,  a zero-crossing 
estimator,  which  has  been  very  popular,  has  been  used  to  estimate  the  mean  frequency. 

The  phase  detector.  Equation  (2-12)  is  an  expression  for  the  mean  frequency  of  a 
continuous  signal.  In  the  phase  detector,  the  derivatives  of  the  quadrature  components  are 
replaced  by  the  backward  difference  form  (f'(t)  «f(t)-f(t-  At)).  The  mean  angular 
frequency  in  Eq.  (2-12)  then  becomes 
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Z (x(k-l)y(k)-x(k)y(k-l)) 

k=l 

Z (xJ(k)  + y!(k» 

k=l 


(2-18) 


An  algorithm  based  on  the  above  equation  is  called  a phase  detector.  This  is  also  referred 
to  as  the  I/Q  algorithm,  because  it  derives  the  mean  frequency  density  from  the  Doppler 
in-phase  and  quadrature  signal.  This  method  has  a frequency  estimation  range  of -0.22  to 
0.25  PRF  and  was  found  to  be  a biased  estimator  [Barber  et  al.,  1985]. 

The  instantaneous  frequency  estimator.  The  second  method  utilizes  the  fact  that  the 
mean  frequency  is  the  time  average  of  the  instantaneous  frequency.  The  numerator  in  Eq. 
(2-12)  can  be  expressed  in  terms  of  the  instantaneous  frequency.  Using  the  definition  of 
the  instantaneous  frequency,  we  have 

x(t)  y'(t)  - x'(t)  y(t)  = e2(t)  co,(t)  (2-19) 


where  e^t)  = x2(t)  + y2(t).  Using  Eq.  (2-19)  and  Eq.  (2-12),  we  obtain 

00 

/ e2(t)  ©,(t)  dt 

com  = ^ (2-20) 

/ e2(t)  dt 


For  the  discrete  signal,  the  mean  angular  frequency  becomes 


^e2(k)©;(k) 
k iyoo 


(2-21) 


The  instantaneous  frequency  in  the  numerator  was  calculated  based  on  Eq.  (2-5).  For  the 
discrete  signal  the  instantaneous  angular  frequency  can  be  written  as 
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co , (k)  - [9(k)  - 0(k  - 1)]  - tan-'  ^ - tan-'  . (2-22) 

x(k)  x(k-l) 

Using  the  identity  for  the  difference  of  the  arc  tangents,  and  rearranging  terms,  Eq.  (2-22) 
becomes 


co  (1c)  = tan' 

x(k)x(k-l)  + y(k)y(k-l) 


(2-23) 


Mean  frequency  estimation  based  on  Eq.  (2-21)  together  with  Eq.  (2-23)  is  referred  to  as 
the  instantaneous  frequency  algorithm,  because  the  mean  frequency  is  derived  from  an 
appropriately  weighted  average  of  instantaneous  frequencies.  Frequency  estimation 
ranges  from  -0.5  to  0.5  PRF.  For  signals  without  amplitude  modulation,  the  instantaneous 
method  produced  unbiased  estimation  for  the  mean  frequency  [Barber  et  al.,  1985], 

The  complex  autocorrelator.  As  shown  in  Eq.  (2-10),  the  mean  frequency  can  be 
expressed  in  terms  of  the  ACF.  Kasai  et  al.  [1985]  proposed  a method  to  estimate  the 
mean  frequency  and  the  variance  of  the  spectrum  using  the  ACF.  As  shown  in  Eq.  (2-10) 
and  (2-17),  the  mean  frequency  and  the  variance  of  a spectrum  can  be  written  in  terms  of 
the  ACF  and  derivatives  of  the  ACF 


,1 

j Rzz(0)  Rxx(0)+Ryy(0) 


(2-10) 


CT2  = 


R'g(Q) 

Rzz(°) 


R,,zz(°) 

RJ0) 


(2-17) 


If  we  express  the  R^x)  as 

RJ*)  =|RZz(T)lexP(j<J)('0)  (2-24) 

and  considering  that  |R2z(x)|  is  an  even  function  and  <[)(x)  an  odd  function,  we  have 
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R,=  (x)  = (|Rzz(t)r+j  |R„(t)|  exp(j<j)'(T)). 


(2-25) 


Therefore,  when  x=0,  we  have 


R'=(0)=j|Rjo)|  *•(<>) 


(2-26) 


rjo)-I*jo)I 


(2-27) 


Inserting  Eq.  (2-26)  and  (2-27)  into  (2-10), 
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= <t>'(0). 


(2-28) 


For  discrete  signals,  we  replace  <J>'(0)  with  <|>(1)  - <j>(0),  therefore  the  mean  angular 
frequency  becomes 


The  variance  is  an  approximation  of  a Taylor  series.  The  implementation  of  this  method 
requires  signal  power  and  a complex  ACF  of  one  sample  lag.  Kasai  et  al.  used  a thread 
phantom  to  evaluate  their  method,  and  found  a linear  relationship.  This  method  was 
applied  successfully  to  real-time  two-dimensional  blood  flow  imaging  where  a fast 
computation  is  mandatory. 

Zero-crossing  method.  The  simplest  approach  in  estimating  the  mean  frequency  of 
the  Doppler  spectrum  is  based  on  the  density  of  the  zero-crossings.  This  method  has  been 
applied  in  most  Doppler  systems  in  the  continuous  and  pulsed  wave  modes.  In  the  case  of 
a continuous  signal,  directional  zero-crossing  detection  was  performed  by  averaging  short 
pulses  generated  at  the  positive  crossing  of  one  of  the  quadrature  signals  at  a threshold 
level,  while  the  sign  of  the  pulse  is  determined  by  the  polarity  of  the  other  quadrature 


= <K°)  = 4>(1)  * 4>(°)  = 4>(1)- 


(2-29) 


Similarly,  the  variance  can  be  written  as 


(2-30) 
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signal.  This  concept  can  be  modified  for  a discrete  signal.  The  output  of  the  time-discrete 
directional  zero-crossing  detector  can  be  expressed  as 

f (k)  = (-sgn(x(k))  + sgn(x(k  - l))(sgn(y(k))  + sgn(y(k  - 1)))  / 8 
+ (-sgn(x(k))  + sgn(x(k  - l))(sgn(y(k))  + sgn(y(k  - l)))sgn(<  f (k)  >)  / 8. 

The  function  sgn(-)  indicates  the  signum  function  being  -1  if  the  argument  is  negative,  zero 
if  the  argument  is  zero,  and  +1  if  the  argument  is  positive.  The  first  term  of  the  right-hand 
side  of  the  above  expression  indicates  the  detected  zero-crossings,  the  polarity  of  which 
could  be  established  on  the  basis  of  the  phase  relationship  between  both  quadrature 
signals.  The  second  term  contains  the  detected  zero-crossing  with  ambiguity  about  the 
polarity. 

The  zero-crossing  method  introduces  a bias  in  the  estimated  mean  frequency,  which 
is  related  to  the  bandwidth  of  the  Doppler  signal  [Van  Leeuwen  et  al.,  1986],  This  error  is 
less  significant  in  a high-resolution  PW  Doppler  system  in  which  the  bandwidth  of  the 
Doppler  signals  is  generally  narrow. 

2.4.2.3.  Performance  of  the  estimators 

The  performance  of  the  previously  mentioned  mean  frequency  estimators,  in  terms  of 
the  bias  and  variance,  is  affected  by  several  factors  including  noise  and  bandwidth.  When 
measuring  the  blood  velocity  deep  inside  the  human  body,  the  returned  ultrasound  echoes 
are  hampered  by  a low  signal-to-noise  ratio  (SNR).  The  spectral  shape  and  the  bandwidth 
of  the  Doppler  signal  also  influence  the  behavior  of  the  estimators.  The  spectral  shape  and 
the  bandwidth  are  mainly  affected  by  the  velocity  profile  and  the  insonating  angle. 
Therefore,  there  could  be  variance  between  various  arteries  and  various  individuals. 
Aliasing  could  be  another  source  of  estimation  error  depending  on  the  implementation 
techniques. 
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Many  experimental  results  have  been  published  on  the  performance  of  the  time- 
domain  Doppler  frequency  estimators.  Most  performance  results  are  described  in  terms  of 
the  bias  and  variance  of  the  estimate  with  respect  to  various  noise  levels  and  bandwidths. 
Hoeks  et  al.  [1984]  studied  the  performance  of  the  instantaneous  and  zero-crossing 
estimators.  Barber  et  al.  [1985]  worked  on  the  phase  and  instantaneous  frequency 
estimator.  Kasai  et  al.  [1985]  proposed  the  autocorrelation  estimator.  Van  Leeuwen  et 
al.  [1986]  tested  the  performance  of  all  four  estimators.  The  researchers  focused  the  use 
of  the  zero-crossing  method. 

There  have  been  several  ways  proposed  to  validate  the  frequency  estimators.  The 
tested  Doppler  signal  is  collected  either  from  a computer  simulated  signal  [Van  Leeuwen 
et  al.,  1985]  or  flow  phantom  [Kasai  et  al.,  1985]  or  real  blood  flow  [Reneman  et  al., 

1973].  When  a Doppler  signal  is  collected  from  real  blood  flow,  other  measurement 
techniques,  such  as  the  electromagnetic  flowmeter,  are  employed  in  order  to  validate  the 
flow  information  from  the  time-domain  frequency  estimation. 

The  phase  detector.  The  phase  detector  is  a biased  estimator  and  has  limited 
frequency  estimation  range  [Barber  et  al.,  1985],  This  problem  can  be  illustrated  by  the 
substitution  of  the  monocomponent  signal  x(k)  = exp(jco0k),  yielding  co^k)  = sin(co0) 
instead  of  the  true  frequency  of  co0.  Van  Leeuwen  et  al.  [1986]  reported  that  only  for 
relatively  narrowband  signals  and  a favorable  SNR,  an  acceptable  error  was  found  for  the 
mean  frequency  up  to  0.2  PRF.  Increasing  the  bandwidth  or  decreasing  the  SNR  resulted 
in  a rapid  deterioration  of  the  estimator  performance. 

The  instantaneous  frequency  estimator.  The  instantaneous  frequency  estimator  is  an 
unbiased  estimator  and  has  the  estimation  frequency  range  of  the  Nyquist  interval,  i.e., 
from  -0.5  PRF  to  0.5  PRF  [Barber  et  al.,  1985].  Hoeks  et  al.  [1984]  provided  the 
theoretical  expression  for  the  probability  distribution  of  the  instantaneous  frequency.  They 
reported  underestimation  of  the  mean  frequency  for  wideband  signals  and  for  signals  with 
a low  SNR.  This  underestimation  can  be  explained  as  follows:  The  instantaneous 
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frequency  estimate  can  have  the  range  wider  than  half  of  the  sampling  rate  (the  Nyquist 
interval).  Therefore,  it  is  likely  that  correctly  sampled  signals  will  generate  instantaneous 
frequencies  outside  the  Nyquist  interval.  However,  due  to  the  sampling  procedure  these 
instantaneous  frequencies  will  be  mapped  within  the  Nyquist  interval,  causing  an 
underestimation  of  the  mean  frequency,  especially  for  the  mean  frequency  close  to  the 
boundaries  of  the  Nyquist  interval. 

Barber  et  al.  [1985]  found  that  the  performance  of  the  instantaneous  frequency 
method  is  superior  except  near  zero  velocity,  where  the  phase  detector  prevails.  For  this 
reason,  they  employed  a hybrid  method  where  the  high  velocity  was  estimated  by  the 
instantaneous  frequency  method  and  low  velocity  was  estimated  by  the  phase  detector. 
They  compared  the  performance  of  the  hybrid  method  with  a spectral  analysis  method 
using  the  FFT.  In  the  spectral  analysis,  the  power  spectral  centroids  were  used  as  a mean 
frequency  estimate.  They  reported  that  in  poor  SNR,  the  performance  of  the  hybrid 
method  was  superior  to  the  spectral  method. 

The  complex  autocorrelator.  The  complex  autocorrelator  requires  the  complex  ACF 
of  an  analytic  signal  at  lag  1,  i.e.,  11^(1).  The  main  advantage  of  the  method  is  that  it  is 
not  hampered  by  the  effect  of  the  averaging  circular  spectrum  problem  which  causes 
underestimation  in  case  of  the  instantaneous  frequency  method.  In  order  to  explain  how 
this  is  possible,  let's  consider  a continuous,  but  partially  aliased,  wideband  PSD  S^co)  in 
Fig.  2-3.  The  corresponding  ACF  of  S^co)  is  defined  as  the  inverse  Fourier  Transform 


(2-32) 


For  m = 1, 
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S(co) 


Fig.  2-3  An  example  of  power  spectral  density 
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(2-33) 


Using  the  fact  that  b = 2k  + c and  8^(0+271)=  S^co),  R^/l)  becomes 


(2-34) 


Therefore,  the  autocorrelation  method  overcomes  the  problem  of  circular  spectrum 
averaging. 

Another  advantage  of  the  autocorrelator  is  that  adding  white  noise  does  not  affect  the 
complex  ACF  11^(1).  If  a signal  is  contaminated  with  an  additive  white  noise  with  power 
of  N0,  and  if  the  signal  and  the  noise  are  independent,  we  can  write 


where  S%(co)  and  SD(co)  are  PSD's  of  measured  and  true  signals,  respectively.  Then,  the 
ACF  R^(l)  becomes 


which  is  the  same  as  R^l).  Therefore,  the  autocorrelator  can  achieve  robust  estimation 
of  the  mean  frequency. 

Unlike  the  previous  frequency  estimator,  the  autocorrelator  provides  an  adequate 
estimate  of  the  mean  frequency  of  the  Doppler  signal,  in  dependent  of  noise  level, 
although  the  SNR  will  affect  the  variance  of  the  estimate  [Van  Leeuwen  et  al.,  1986]. 
Kasai  et  al.  [1985]  reported  that  frequency  estimation  of  the  Doppler  signal  from  a string 
phantom  provided  a linear  relationship  with  the  string  velocity.  Using  a simulation  signal, 


S^(co)  = S»  + N0, 


(2-35) 


(2-36) 
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Van  Leeuwen  et  al.  [1986]  concluded  that  the  autocorrelator  yields  the  best  performance 
among  the  three  other  estimators. 

Zero-crossing  method.  Flax  et  al.  [1970]  investigated  the  performance  of  the  zero- 
crossing for  a CW  ultrasound  Doppler.  They  showed  that,  in  the  presence  of  a broad 
spectrum,  the  variance  of  the  flow  estimate  increased.  Since  the  Doppler  spectrum  shape 
is  dependent  of  the  velocity  profile,  the  performance  of  the  zero-crossing  is  dependent  on 
the  blood  velocity  profile.  For  instance,  Woodcock  et  al.  [1972]  reported  that,  for  a 
parabolic  profile,  the  zero-crossing  read  16%  higher  that  the  mean  velocity;  for  a flat 
profile,  the  reading  was  an  accurate  mean  value.  This  implies  that  there  is  a tremendous 
problem  in  estimating  the  mean  blood  flow  at  an  arbitrary  artery,  since  the  blood  velocity 
profile  differs  among  the  various  arteries  and  during  a heart  beat.  In  a clinical  vascular 
laboratory  study,  Jonhston  et  al.  [1977]  found  that  a zero-crossing  detector  system 
accurately  recorded  the  velocity  wave  in  only  15%  of  the  femoral,  3 1%  of  the  popliteal, 
67%  of  the  dorsalis  pedis,  and  62%  of  the  posterior  tibial  artery  recordings. 

The  ultrasound  Doppler  signal  is  affected,  not  only  by  the  movement  of  RBCs,  but 
also  by  the  movement  of  vessel  walls.  Wall  velocities  are  usually  lower  than  RBC 
velocities,  and  therefore  are  associated  with  Doppler  frequencies  less  than  100  Hz  [Flax  et 
al.,  1973],  In  addition,  the  Doppler  signal  by  the  wall  motion  is  approximately  1-2  orders 
of  magnitude  larger  than  the  flow-induced  signal.  Therefore,  wall-motion  signals  in  effect 
lift  the  flow-induced  signal  above  the  zero  axis  so  that  many  zero-crossings  are  not 
counted  by  the  zero-crossing  detector. 

2,4.3.  Spectral  Analysis 

The  implementation  of  the  time-domain  analysis  provides  simple  algorithm  so  that  it 
can  be  achieve  faster  processing  time  than  spectral  analysis.  However,  experimental 
results  indicate  that  it  is  inappropriate  for  universal  applications.  Spectral  analysis  involves 
the  estimation  of  short-time  power  spectral  density  (PSD)  of  the  Doppler  signal.  In  the 
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70's,  a filterbank  was  used  in  the  analysis.  [Light,  1970;  Johnston  and  Taraschuk,  1976], 
However,  the  fast  Fourier  transform  (FFT)  algorithm  [Cooley  and  Tukey,  1965]  and  the 
increased  digital  computing  power  allowed  this  technique  to  be  implemented  inexpensive 
and  in  real-time.  The  main  advantage  of  spectral  analysis  is  that  it  can  present  a complete 
picture  of  the  spectral  distribution,  while  the  time-domain  analysis  is  limited  to  the  mean 
and  the  variance.  Spectral  analysis  can  provide  a good  display  of  the  frequency  changes 
during  a cardiac  cycle,  particularly  when  the  signal-to-noise  ratio  is  poor.  Also,  the 
spectral  analysis  allows  us  to  recognize  overlapping  signals  that  may  arise  from  a 
neighboring  vessel.  Spectral  analysis  is  very  robust  in  detecting  turbulence  flow,  because 
turbulent  flow  introduces  a broad  spectrum.  Except  for  certain  locations,  the  turbulent 
blood  flow  indicates  arterial  stenosis. 

Alternative  spectral  estimation  techniques  have  been  applied  to  the  Doppler  signal.  In 
his  series  of  studies  [Kaluzynski,  1987]  Kaluzynski  investigated  the  application  of 
autoregressive  spectral  estimation  (ARSE)  to  the  Doppler  signal  by  comparing  the 
periodogram  and  ARSE  from  steady  flow,  aortic  flow  and  branchial  artery  flow.  He 
reported  that  the  periodogram  yields  a statistically  unstable  estimate  and  does  not  provide 
a reliable  spectrum  even  in  the  case  of  a steady  model  flow.  He  concluded  that  the  AR 
method  gives  the  best  results  provided  the  spectrum  is  narrow,  even  for  short  data  lengths. 
In  his  second  study  [Kaluzynski,  1989]  he  applied  periodogram,  the  autocorrelation 
method,  the  maximum  likelihood  method  and  the  Wigner-Ville  distribution  to  the  Doppler 
signal  obtained  from  a fully  insonated  laminar  model  flow.  Although  he  concluded  that 
the  maximum  likelihood  method  outperformed  the  others,  the  experimental  setting  of  the 
steady  laminar  model  flow  is  fat  removed  from  the  conditions  in  human  body.  Therefore, 
application  of  his  results  to  the  human  artery  would  impose  some  limitations. 

Vaitkus  and  Cobbold.  [1988a,  1988b]  investigated  the  performance  of  the  available 
parametric  spectral  estimation  methods  using  a simulated  Doppler  signal  that  was  derived 
from  a spectral  shape  using  the  inverse  Fourier  transform.  The  simulated  spectral  data 
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was  previously  shown  to  be  very  close  to  that  obtained  by  ensemble  averaging  of  a normal 
carotid  waveform  at  around  peak  systole  [Mo  and  Cobbold,  1986a],  After  the  spectral 
estimation,  they  determined  the  performance  of  each  estimation  method  using  1)  bias, 
variance  and  mean  square  error,  2)  two  spectral  distance  indices,  and  3)  spectral  flatness 
index.  They  found  that  AR  (p=8,  autocorrelation  method)  and  ARMA  (p=4,  q=4,  singular 
value  decomposition)  showed  good  agreement  with  the  theoretically  simulated  spectrum. 

2,4.4.  Velocity  Envelop  Processing 

Even  though  the  Doppler  spectra  contain  the  complete  information  of  the  blood 
velocity,  most  clinical  applications  of  the  UDS  are  interested  only  in  a part  of  the  Doppler 
spectra.  In  most  clinical  applications,  the  UDS  is  used  as  a flowmeter.  For  example,  we 
can  use  a UDS  to  monitor  the  mean  blood  flow  to  the  human  brain  during  general 
anesthesia.  Many  clinically  significant  features  can  be  derived  from  the  instantaneous 
blood  flow  information.  Researchers  have  investigated  various  techniques  used  to  extract 
blood  flow  information  from  the  Doppler  spectra.  The  most  frequently  used  method 
estimates  the  maximum  frequency,  fj^.  In  this  method,  the  following  conditions  are 
assumed. 

1)  The  true  Doppler  signal  is  low-pass  stochastic  process 

S(f)  = 0,forf>fm4X  (2-37) 

where  S(f)  is  the  signal  power  spectral  density  (PSD)  of  s(n). 

2)  The  observed  Doppler  signal,  x(n),  is  the  sum  of  the  true  Doppler  signal,  s(n),  and 
white  Gaussian  noise,  N(n),  that  is, 

x(n)  = s(n)  + N(n).  (2-38) 

3)  The  noise  is  not  correlated  with  the  true  Doppler  signal,  which  allow  us  to  express 
the  PSD  of  x(n),  X(f),  such  that 
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X(f)  = S(f)  + N0  (2-39) 

where  N0is  the  noise  PSD. 

To  date,  most  maximum  frequency  estimators  are  based  on  spectral  analysis  by  FFT. 
Insights  on  the  various  maximum  frequency  estimators  can  be  gained  by  defining  the 
integrated  Doppler  PSD  such  that 


4>(f)=£x(k) 

k=fL 


(2-40) 


where  X(k)  is  calculated  by  FFT. 

Mo  et  al.  [1988]  investigated  the  performance  of  four  methods  for  estimating  the 
maximum  frequency  envelope  from  the  ultrasound  Doppler  spectrum.  They  are: 
percentile  method,  D'Alessio's  threshold  crossing  method,  a modified  threshold  crossing 
method,  and  a new  hybrid  method. 

In  the  percentile  method,  the  maximum  frequency  estimate  is  defined  as  f^  such  that 


100 -a 
100 


<J)T 


(2-41) 


where  <{>T  is  the  total  signal  plus  noise  power,  and  a is  a chosen  parameter  whose  optimum 
value  depends  on  both  the  SNR  and  the  signal  bandwidth. 

D'Alessio  proposed  an  objective  algorithm  for  the  maximum  velocity  estimation  based 
on  the  statistical  characteristics  of  the  FFT  spectral  estimator.  He  assumed  that  the 
observed  signal  is  the  summation  of  the  Doppler  signal  and  noise.  Furthermore,  he 
assumed  that  the  Doppler  signal  was  modeled  as  a stationary  zero-mean  Gaussian  process 
and  that  the  noise  was  a white  Gaussian  process.  Starting  from  the  noise  end  of  the  total 
spectrum,  spectral  samples  in  successive  frequency  bins  are  compared  to  a threshold  level. 
If,  in  a sequence  of  r successive  bins,  there  are  at  least  m bins  that  exceed  the  threshold, 
X^,  then  the  maximum  frequency  estimate  is  taken  to  be  the  highest  bin  frequency  in  that 
sequence.  For  example,  when  r = m = 2,  the  threshold  level  is  given  by 
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X^  = -N0ln(l-P,)  (2-42) 

where  Ps  is  a chosen  minimum  specificity  setting. 

D'Alessio's  algorithm  is  based  on  the  Doppler  spectrum  using  a rectangular  window. 
Mo  et  al.  argued  that  a similar  threshold  method  can  be  used  with  a nonrectangular 
window  even  though  the  statistic  of  the  noise  spectrum  would  be  altered  so  that 
D'Alessio's  algorithm  could  no  longer  be  applied.  This  method  was  called  as  the  modified 
threshold  crossing  method. 

Finally,  Mo  et  al.  proposed  a hybrid  method  combining  the  features  of  both  the 
percentile  method  and  the  threshold  method.  The  maximum  frequency  fj^  is  given  by  the 
lowest  frequency  in  the  interval  [fL,  fH]  that  satisfies 

<Kfmax  ) = <i>T  - X™  (fH  ~ fm«  ) (2-43) 

where  X”{  is  chosen  such  that  it  is  greater  than  or  equal  to  N0 

For  the  threshold  and  the  hybrid  methods,  the  noise  power  should  be  estimated.  In 
practice,  the  noise  power  is  estimated  by  averaging  over  the  tail  end  of  the  PSD  of 
measured  signal.  The  rationale  of  the  percentile  method  is  that  the  power  of  the  white 
noise  is  the  same  during  the  measurement,  and  the  constant  a implies  the  magnitude  of  the 
noise  power.  The  threshold  crossing  and  hybrid  methods  require  to  estimate  the  noise 
power  so  that  the  maximum  velocity  is  adaptive  to  the  noise. 

The  performances  of  the  four  maximum  frequency  estimators  have  been  assessed  and 
compared  using  a simulated  signal.  Mo  et  al.  reported  that  the  modified  threshold 
crossing  and  the  hybrid  methods  have  the  best  and  most  consistent  performance  over  the 
wide  range  of  conditions  tested. 

2.5.  Medical  Applications-Feature  Extraction  and  Hemodynamics 
The  practical  objective  of  the  signal  processing  of  the  ultrasound  Doppler  is  to 
develop  and  identify  an  efficient  method  to  extract  features  of  interest.  The  Doppler  signal 
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is  further  processed  filtering  out  inappropriate  or  irrelevant  information.  There  is  no 
agreement  as  to  the  best  method  for  extracting  these  features.  A rule  of  thumb  is  to 
choose  the  method  based  on  the  intended  application  of  the  ultrasound  Doppler. 

In  this  section,  we  first  review  factors  that  affect  blood  flow.  This  is  the  main 
objective  of  hemodynamics.  Hemodynamics  studies  the  forces  generated  by  the  heart  and 
the  resulting  motion  of  the  blood  through  the  cardiovascular  system  [Milnor,  1989], 
Research  in  hemodynamics  is  further  divided  into  the  following  three  areas 

1)  the  physical  properties  of  the  heart,  blood,  and  blood  vessels; 

2)  the  relation  of  these  properties  to  the  phenomena  that  are  observed  in  the 
circulation; 

3)  and  the  application  of  the  results  to  physiological  research  or  clinical  medicine. 

The  ultrasound  Doppler  can  be  considered  as  a tool  to  observe  blood  circulation.  In  the 
application  of  the  ultrasound  Doppler,  the  following  hemodynamic  factors  are  of  particular 
interest:  amount  of  blood  flow  and  control  mechanism  of  the  blood  flow;  pulsatile  nature 
of  the  arterial  blood  flow;  and  flow  disturbances  associated  with  vascular  disease  such  as 
stenosis.  We  will  discuss  how  the  hemodynamic  factors  are  mapped  into  the  Doppler 
signal,  and  how  to  design  the  feature  extraction  phase  in  order  to  retrieve  the 
hemodynamic  parameters.  Next,  we  will  consider  other  ultrasound  instrumentation  that 
can  provide  hemodynamic  information. 

2.5.1,  Hemodynamics  and  Doppler  Signal 

The  blood  supplies  nutrient  and  oxygen  throughout  the  human  body.  The  amount  of 
blood  supply  to  an  organ  depends  on  two  factors:  perfusion  pressure  and  the  vascular 
resistance.  The  perfusion  pressure  is  the  pressure  difference  between  the  arterial  pressure 
and  the  terminal  blood  pressure  of  the  organ.  The  vascular  resistance  is  a lumped 
parameter  which  is  the  function  of  the  viscosity  of  the  blood,  the  total  vessel  diameter,  and 
the  length  of  the  blood  vessel.  Figure  2-4  depicts  an  analogy  of  the  human  circulation  to 
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Resistance 


Fig.  2-4  An  analogy  of  blood  circulation  to  an  electric  circuit 
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an  electrical  circuit.  The  human  body  has  a complex  control  mechanism  to  distribute  the 
blood  to  each  organ.  One  way  to  control  the  blood  flow  is  to  change  the  vascular 
resistance.  When  a part  of  the  human  body  needs  blood,  blood  vessels  supplying  blood 
throughout  that  part  dilate.  As  a result,  the  resistance  decreases  and  the  blood  flow 
increases.  This  is  called  vasodilation.  A reverse  action  is  executed  when  it  wants  to 
reduce  the  blood  supply.  One  of  the  main  indicators  of  the  resistance  of  the  distal  vascular 
bed  is  the  diastolic  blood  flow.  For  example,  the  internal  carotid  artery,  because  it 
supplies  low-resistance  capillaries,  has  larger  diastolic  flow  than  the  external  carotid  artery 
with  high  resistance  capillaries. 

One  of  the  problems  of  the  UDS  is  that  it  measures  the  blood  velocity  rather  than 
blood  flow.  In  order  to  measure  the  blood  flow  we  have  to  know  the  diameter  of  the 
blood  vessel.  Additional  instrumentation  is  needed  to  measure  the  dimensions  of  the 
blood  vessels.  Nevertheless,  if  we  assume  that  the  diameter  of  the  blood  vessel  is  constant 
during  a cardiac  cycle,  the  ultrasound  Doppler  can  be  used  as  a blood  flowmeter.  The 
assumption  is  valid  when  the  blood  velocity  is  measured  at  a conducting  artery.  The 
conducting  arteries  have  thick  vessel  walls  so  that  we  can  assume  a constant  vessel 
diameter. 

One  of  the  most  important  properties  of  blood  flow  is  that  it  is  pulsatile.  The  most 
prominent  studies  on  the  pulsatile  blood  flow  are  by  Womersley  and  McDonald 
[McDonald,  1974],  Womersley  derived  the  shape  of  the  blood  flow  distribution.  Milnor 
[1989]  reviewed  the  various  assumptions  that  were  made  in  the  Womersley  analysis.  A 
typical  behavior  of  the  flow  distribution  can  be  observed  within  a cardiac  cycle.  During  an 
acceleration,  the  flow  distribution  tends  to  be  flat  with  all  of  the  fluid  acceleration  at 
approximately  the  same  velocity.  After  the  peak  velocity  is  reached,  the  center  stream 
tends  to  speed  ahead  of  the  fluid  near  the  walls  because  the  flow  near  the  walls  is  slowed 
by  viscous  friction. 
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The  Doppler  signal  is  influenced  by  the  pulsatile  nature  of  the  blood  flow.  Since  the 
Doppler  signal  is  modeled  as  a stochastic  process,  the  pulsatile  blood  flow  results  in  a 
nonstationary  Doppler  signal.  The  distribution  of  blood  flow  in  a sample  volume 
determines  the  bandwidth  of  the  spectrum  of  the  Doppler  signal.  Two  more  effects  can 
influence  the  bandwidth  irrespective  of  the  blood  flow  distribution.  One  is  the  transit  time 
effect  of  the  red  blood  cells  moving  through  the  ultrasound  beam.  The  transit  time  is  the 
interval  during  which  each  red  blood  cell  stays  in  a sample  volume.  As  the  blood  velocity 
is  increased,  the  transit  time  decreases  and  the  bandwidth  of  the  Doppler  spectrum 
increases.  The  other  is  the  geometrical  broadening  due  to  the  diameter  of  the  transducer 
in  relation  to  the  wavelength  of  the  transmitted  ultrasound.  The  effect  strongly  depends 
on  the  distance  and  will  be  manifested  when  the  probe  is  close  to  the  sampling  volume. 

In  laminar  flow,  the  movement  of  a single  cell  coincides  with  an  axis.  However,  in 
turbulent  flow,  the  movement  of  blood  cell  is  random.  The  transition  from  laminar  to 
turbulent  flows  occurs  when  the  Reynolds  number  exceeds  a threshold  value  [Milnor, 
1989].  Also,  in  the  presence  of  stenosis,  the  poststenotic  flow  becomes  turbulent.  In  the 
normal  human  artery,  blood  flow  is  laminar  with  a few  exceptions.  The  turbulence  flow 
due  to  the  presence  of  stenosis  tends  to  broaden  the  Doppler  spectrum  during  a systolic 
period. 

2.5.2,  Feature  Extraction 

Feature  extraction  is  based  on  a waveform  derived  by  time-domain  analysis  or 
spectral  analysis  followed  by  a velocity  envelope  estimator.  This  waveform  contains  the 
information  of  blood  flow.  In  the  later  case,  there  are  three  ways  to  calculate  the  velocity 
envelope  from  the  Doppler  spectrum:  maximum,  mean,  and  mode  velocity  methods. 
However,  maximum  velocity  envelope  method  has  been  shown  to  be  the  most  robust  and 
sensitive  way  to  extract  quantitative  features  [Kirkham  et  al.,  1986;  Thompson  et  al., 
1986b],  The  first  step  in  feature  extraction  is  to  identify  each  cardiac  cycle.  This  can  be 
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done  by  detecting  the  minimum  and  the  maximum  of  the  waveform  or  using  other  signals 
such  as  the  electrocardiogram  (ECG).  Next,  three  basic  values  are  calculated:  the 
maximum,  the  minimum,  and  the  mean  value.  In  the  medical  literature,  the  maximum  and 
minimum  are  called  the  systolic  and  diastolic  velocities,  respectively. 

The  mean  of  the  waveform  is  used  as  a direct  indicator  for  the  blood  flow.  In  a 
clinical  experiment  using  a Transcranial  Doppler  (TCD)  [Kirkham  et  al.,  1986],  the  mean 
of  the  maximum  velocity  envelope  was  calculated  from  the  sample  volume  at  the  Ml 
section  of  a middle  cerebral  artery  (MCA)  and  at  the  internal  carotid  artery  (ICA).  The 
relationship  between  end-expiratory  pC02  and  mean  of  maximum  envelope  was  found  to 
be  linear.  Furthermore,  the  slope  of  the  linearity  was  very  similar  to  those  reported  using 
direct  methods  of  measuring  cerebral  blood  flow.  This  fact  may  imply  that  it  is  unlikely 
that  a change  in  diameter  of  either  the  ICA  or  the  MCA  occurred  when  the  pC02  level  is 
changing.  Therefore,  the  blood  velocity  may  be  used  to  monitor  changes  in  flow. 

There  are  many  other  quantitative  features  that  characterize  the  velocity  envelope  and 
have  a relation  with  hemodynamic  properties.  We  will  review  the  pulsatility  index  (PI), 
the  resistance  index  (RI),  the  and  spectral  broadening  index  (SBI). 

The  pulsatility  index  (PI)  is  defined  as 


VE  - VF 

pj  _ YJ^max  V JJ/min 

VE 


(2-44) 


where  VE^  denotes  the  maximum  value  of  the  velocity  envelope,  and  so  on.  The  PI  is 
known  to  be  useful  for  detecting  occlusive  arterial  disease  [Johnston  et  al.,  1978;  Evans  et 
al.,  1981;  Johnston  et  al.,  1983;  Thompson  et  al.,  1986a  and  1986b;  Thompson  and 
Trudinger,  1990],  The  resistance  index  (RI)  is  defined  as 
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(2-45) 


The  RI  is  also  called  the  Pourcelot  ratio.  This  has  been  used  to  determine  vascular 
impedance  [Fitzgerald  et  al.,  1984;  Thompson  et  al.,  1986a  and  1986b]  Note  that  the  PI 
and  RI  are  independent  of  the  angle  between  the  probe  and  the  vessel,  since  the  unknown 
constant  of  proportionality  is  present  in  both  the  numerator  and  the  denominator  and 
disappears  in  the  calculation.  The  PI  and  RI  are  defined  by  a velocity  envelope  waveforms 
within  a cardiac  cycle.  The  spectral  broadening  index  (SBI)  is  defined  as 


where  and  fmem  are  the  maximum  and  mean  frequency  of  the  Doppler  spectrum  at  a 
given  time  [Wijn  et  al.,  1987],  SBI  has  been  used  especially  during  systole  but  for  each 
time  SBI  can  be  derived  from  the  max  and  mean  curves. 


There  are  five  noise  sources  of  interest. 

1)  weak  ultrasound  return  from  a sample  volume 

2)  interference  by  a subject 

3)  movement  of  the  blood  vessel 

4)  interference  from  other  devices 

5)  errors  from  the  Doppler  signal  processing 

The  first  is  caused  by  a weak  scattered  ultrasound  wave  from  a sample  volume.  This 
can  be  caused  by  high  absorption  or  reflections  by  the  media  between  the  probe  and  the 
sample  volume.  It  has  been  known  that  ultrasound  does  not  penetrate  human  scull  well 
with  higher  age  group.  In  some  cases,  Doppler  signal  can  not  be  detected  regardless  of 
sex  and  age.  The  strength  of  the  signal  can  be  represent  by  a signal-to-noise  ratio  (SNR). 


oJol  = 

f 

max 


CDT max  mean 


(2-46) 


2.6,  Noise  Analysis 
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The  second  type  of  noise  is  interference  from  other  devices.  Operations  in  a noisy 
environment  such  as  the  operating  room  can  deteriorate  the  quality  of  the  Doppler  signal. 
The  third  noise  arises  from  the  subject  to  whom  the  probe  is  attached.  During  the 
intraoperative  monitoring,  head  gear  is  mounted  on  a patient  for  continuous  attachment  of 
the  probe.  However,  voice  and  head  movement  by  the  patient  contaminate  the  Doppler 
signal.  The  fourth  is  due  to  vibration  of  the  blood  vessel.  Most  TCD  provide  a highpass 
filter  to  eliminate  the  assumed  low  frequency  noise. 

Noise  types  can  be  classified  according  to  their  spectral  contents  and  duration.  The 
first  noise  can  be  considered  a white  long-duration  noise.  Interference  by  voice  appears  as 
the  trace  of  formant  frequencies.  The  noise  from  other  devices  can  be  a continuous 
narrowband  or  intermittent  short-duration  white  noise.  Even  though  the  magnitude  of  the 
blood  vessel  vibrations  is  substantial,  the  vibration  is  largely  within  100  Hz. 

2.1.  Diagnostic  Instruments  Using  Ultrasound 

In  addition  to  the  ability  to  detect  blood  velocity,  ultrasound  can  be  used  to  provide 
other  information.  In  this  section,  we  will  briefly  review  the  available  diagnostic 
instruments  using  ultrasound. 

Although  the  Doppler  spectrum  depends  on  the  velocity  distribution  in  a sample 
volume,  it  is  impossible  to  map  the  Doppler  spectrum  into  a two  dimensional  velocity 
distribution.  However,  a PW  Doppler  with  a smaller  sample  volume  and  the  capability  of 
multigating  can  map  the  velocity  profile  [Histand  et  al.,  1973],  In  order  to  measure  the 
dimension  of  an  artery,  B-mode  ultrasonography  is  used.  Using  the  reflection  properties 
of  ultrasound  between  two  different  media,  it  can  produce  a two  dimensional  image  of  the 
tissue.  Duplex  sonography  combines  the  B-mode  ultrasonography  and  ultrasound 
Doppler.  This  can  produce  the  velocity  information  at  a specific  area  indicated  in  the 
image  of  tissue.  Ultrasound  duplex  color  arteriography  produces  a real-time  color- 
mapped  velocity  information  overlapping  the  two-dimensional  image  of  an  artery.  This  is 
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the  most  recent  and  potentially  one  of  the  most  significant  implementations.  Medical 
applications  of  the  ultrasound  duplex  color  arteriography  are  discussed  in  Merritt  [1990]. 

Transcranial  Doppler  (TCD)  has  been  used  to  detect  the  blood  flow  inside  the  major 
arteries  of  the  human  brain  and  to  measure  the  blood  velocity.  Aaslid  et  al.  [1982]  first 
reported  the  use  of  the  TCD  in  humans,  using  a Doppler  probe  with  a frequency  of  1-2 
MHz  (versus  the  5-10  MHz  of  other  clinically  used  Dopplers)  to  decrease  the  attenuation 
caused  by  bone  and  soft  tissue.  The  axial  length  of  the  sample  volume  is  about  7-13  mm 
which  is  relatively  long  compare  to  other  ultrasound  Doppler  devices. 

2,8,  Summary 

Research  areas  using  ultrasound  Doppler  systems  (UDS)  were  reviewed  being  lumped 
into  three  categories,  namely,  Doppler  signal  modeling,  Doppler  signal  processing,  and 
medical  applications  of  UDS.  In  Doppler  signal  modeling,  the  main  objective  was  to 
analyze  the  interaction  between  the  ultrasound  and  the  blood.  In  Doppler  signal 
processing,  three  distinctive  ways  to  extract  blood  velocity  were  audio  analysis,  time- 
domain  analysis,  and  spectral  analysis.  Time-domain  analysis  provides  a simpler  and  faster 
implementation  than  the  spectral  domain  technique,  however,  the  information  is  limited  to 
the  mean  and  the  variance  of  the  Doppler  spectrum.  The  spectral  analysis  provides 
additional  details  of  the  spectral  contents  of  the  Doppler  signal.  In  medical  applications, 
we  considered  clinically  important  hemodynamic  factors.  Definitions  of  quantitative 
features  and  their  relationship  with  the  hemodynamic  factors  were  discussed.  Finally,  a 
brief  review  of  diagnostic  devices  using  ultrasound  were  presented. 


CHAPTER  3 


PARAMETRIC  SPECTRAL  ESTIMATION-THEORY 
3.1.  Introduction 

In  this  chapter,  we  describe  the  general  theory  of  parametric  spectral  estimation, 
especially  autoregressive  (AR)  spectral  estimation.  Parametric  spectral  estimation 
involves  a technique  to  estimate  the  power  spectral  density  (PSD)  of  a stochastic  process 
based  on  a linear  rational  model.  The  stochastic  process  is  considered  to  be  the  output  of 
the  linear  rational  model  excited  by  white  noise  input.  Kay  [1988,  p.  106]  described  the 
rationale  of  the  parametric  spectral  estimation 

The  classical  methods  of  spectral  estimation  . . . used  Fourier  transform 
operations  on  either  windowed  data  or  windowed  autocorrelation  function 
(ACF)  estimates.  Windowing  of  data  or  ACF  values  makes  the  implicit 
assumption  that  the  unobserved  data  or  ACF  values  outside  the  window  are  zero, 
which  is  normally  an  unrealistic  assumption.  A smeared  spectral  estimate  is  a 
consequence  of  the  windowing.  Often,  we  have  more  knowledge  about  the 
process  from  which  the  data  samples  are  taken,  or  at  least  we  are  able  to  make  a 
more  seasonable  assumption  other  than  to  assume  the  data  or  ACF  values  are 
zero  outside  the  window.  Use  of  a priori  information  (or  assumptions)  may 
permit  the  selection  of  an  exact  model  for  the  process  that  generated  the  data 
samples,  or  at  least  a model  that  is  a good  approximation  to  the  actual  underlying 
process.  It  is  then  usually  possible  to  obtain  a better  spectral  estimate  by  basing 
it  on  the  model  and  estimating  the  parameters  of  the  model  from  the 
observations. 


A key  point  is  that  the  structure  of  the  model  should  accurately  represent  the  signal 
generating  process.  Therefore,  more  knowledge  of  the  signal  is  always  required  for 
parametric  spectral  analysis.  The  parametric  modeling  of  pulse  wave  (PW)  ultrasound 
Doppler  is  presented  in  chapter  4. 
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The  parametric  spectral  estimation  is  accomplished  by  a three-step  procedure.  The 
first  step  is  to  select  a model.  Three  types  of  models  will  be  introduced  to  accomplish  this 
step.  The  second  step  is  to  estimate  the  parameters  of  the  assumed  model  using  the 
available  data  samples.  The  third  step  is  to  obtain  a spectral  estimate  by  substituting  the 
estimated  model  parameters  in  the  theoretical  PSD  implied  by  the  model.  There  are  many 
ways  to  formulate  the  AR  spectral  estimation,  nonetheless,  they  converge  to  similar 
implementations.  We  will  focus  on  the  formulation  of  linear  prediction.  Other 
formulations  will  be  discussed  briefly.  In  the  final  section,  we  will  consider  other  design 
factors  which  may  affect  the  performance  of  the  AR  spectral  estimation. 

3,2.  Definition  of  Linear  Rational  Models 
A linear  rational  model  has  an  input-output  relationship  described  by  a linear 
difference  equation 

x(n)  = Xa(k)x(n-k)  + G^]b(k)u(n-k),  b(0)  = 1 (3-1) 

k=l  k=0 

where  a(k),  b(k)  and  the  gain  G are  the  parameters  of  the  model,  and  u(n)  and  x(n)  are  the 
input  and  output  signals  of  the  model,  respectively.  This  equation  describes  the  output 
x(n)  in  terms  of  a linear  combination  of  past  outputs  and  present  and  past  inputs.  This 
model  is  also  known  as  the  pole-zero  model  of  the  process.  The  transfer  function  of  this 
model  is  given  by 
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1 + 2 b(k)  zk 

H(z)  = G— if (3-2) 

1 + 2 a(k)  zk 

k=l 

where  X(z)  and  U(z)  are  z-transforms  of  x(n)  and  u(n),  respectively.  The  roots  of  the 
numerator  and  denominator  polynomials  are  called  the  zeros  and  the  poles  of  the  model, 
respectively. 

There  are  two  special  cases  of  the  model  that  are  of  interest 

1)  all-zero  model:  a(k)  = 0, 

2)  all-pole  model:  b(k)  = 0. 

The  all-zero  model  is  known  as  the  moving  average  (MA)  model,  and  the  all-pole  model  is 
known  as  the  autoregressive  (AR)  model  [Makhoul,  1975].  The  pole-zero  model  is 
known  as  the  autoregressive  moving  average  (ARMA)  model.  The  most  popular  of  these 
models  is  the  AR  model.  This  is  because  accurate  estimates  of  the  AR  parameters  can  be 
found  by  solving  a set  of  linear  equations.  For  accurate  estimation  of  parameters  of  the 
ARMA  models,  we  need  to  solve  a set  of  highly  nonlinear  equations.  In  the  following 
discussion,  we  will  focus  on  the  all-pole,  or  AR,  model. 

3,3,  Autoregressive  Parameter  Estimation 


3,3.1.  Introduction 

There  are  several  formulations  to  estimate  the  parameters  of  the  AR  model.  Each 
may  be  based  on  different  assumptions.  The  following  formulations  will  be  considered 

- linear  prediction:  Wiener  filtering;  inner  product  formulation 

- AR  parameter  estimation 

- maximum  likelihood  method 

- maximum  entropy  method 
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- correlation  matching 

- Prony's  method 

- spectral  flatness  formulation. 

As  will  be  seen  later,  all  of  these  formulations  result  in  similar  implementations.  In  this 
chapter,  we  will  examine  in  detail  only  the  linear  prediction  formulation  and  briefly  review 
the  main  ideas  for  the  other  formulations.  We  will  discuss  the  assumptions  and  rationale, 
and  the  relationship  between  them.  Some  formulations  listed  above  require  a theoretical 
value  which  can  be  only  obtained  with  an  infinite  observation  window.  Thus,  we  will 
discuss  problems  with  implementation  when  we  have  a finite  set  of  data.  The  following 
implementation  method  will  be  discussed 

- autocorrelation  method;  Yule-Walker  method;  Levinson-Durbin  algorithm 

- covariance  method 

- lattice  method 

- modified  covariance  method;  forward -backward  approach  [Nuttall,  1976];  least 

squares  method  [Ulrych  and  Clayton,  1976] 

- Burg's  method 

We  will  also  show  how  other  formulations  leads  to  the  similar  implementations. 

3,3.2.  Linear  Prediction 

The  term  linear  prediction  was  first  used  by  Wiener  in  1949.  In  speech  processing, 
linear  prediction  has  become  a universally  accepted  term  referring  to  those  methods  of 
speech  analyses  that  result  in  the  solution  of  a predictor  coefficient  based  upon  solving  a 
set  of  linear  simultaneous  equations.  The  problem  of  the  linear  prediction  is  to  predict  the 
unobserved  sample  x based  on  the  observed  data  set  (x(n-l),  x(n-2),  ...  , x(n-p)},  or 
previous  p samples 


50 


x(“)  = -£a(k)x(n-k).  (3-3) 

k=l 

In  linear  prediction,  the  coefficients  -a(k),  i = 1,  2, ...  ,p  , are  called  the  predictor 
coefficients.  The  prediction  error,  e(n),  between  the  actual  value  x(n)  and  the  predicted 
value  k(n)  is  defined  as 


e(n)  = x(n)-x(n)  = x(n)  + J]a(k)x(n-k).  (3-4) 

k=l 

Since  e(n)  represents  the  error  between  the  actual  sample  and  the  predicted  one,  it  would 
seem  reasonable  to  choose  the  predictor  coefficients  so  that  e(n)  is  minimized  in  some 
manner.  Now,  the  problem  of  predicting  the  next  sample  is  transformed  into  the  problem 
of  estimating  predictor  coefficients  that  minimize  the  prediction  error  in  some  manner. 
With  the  method  of  least  squares,  the  prediction  parameters  are  obtained  as  a result  of  the 
minimization  of  total  squared  error  with  respect  to  each  predictor  coefficient.  The  main 
reason  for  this  choice  of  optimization  criterion  is  that  the  resulting  equations  are  linear, 
and  tractable. 

The  total  squared  error  E is  defined  by 

E = 2V(n)  = £(x(n)-x(n))2  (3.5) 

n=a  n=a 

where  a and  b are  the  index  limits  over  which  error  minimization  occurs.  Inserting  x as  in 
Eq.  (3-3),  and  further  manipulating  the  above  equation. 


(3-6) 


where 
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c(i,  j)  = Z x(n  - 0 x(n  ' j)  • (3‘7) 

n=» 

Eq.  (3-6)  shows  that  the  total  squared  error  E is  in  quadratic  form.  Minimization  of  E is 
obtained  by  setting  the  partial  derivative  of  E with  respect  to  a(k),  k = 1,  2, ...  , p to  zero. 

5E  p 

— — = 2^a(i)c(i,k)  = °,  for  k=  1,2,..., p (3-8) 

oa(k)  i=0 

or  since  a(0)  = 1, 

Za(i)c(i>k) = ~c(0,k),  for k = l,2,...,p.  (3-9) 

i=l 

The  p unknown  predictor  coefficients  a(k)  are  obtained  by  solving  this  set  of  p linear 
simultaneous  equations.  The  known  parameters  c(i,k),  for  i = 0,  2, ...  , p,  and  k = 0,  2, ... 

, p are  defined  from  the  observed  data.  After  we  estimate  the  predictor  coefficients,  the 
minimum  total  error  can  be  found  by  inserting  Eq.  (3-9)  into  Eq.  (3-6), 

Emm  =c(0,0)  + £a(k)c(0,k).  (3-10) 

k=l 

Eq.  (3-9)  and  Eq.  (3-10)  are  called  the  Wiener-Hopf  equations. 

We  left  unspecified  the  range  of  the  summation  interval  [a,  b].  Assuming  that  we 
have  N observed  data  (s(n)}  = (s(0),  s(l), ...,  s(N-l)},  the  autocorrelation  method  is 
defined  by  setting  a = -oo  and  b = oo,  and  defining  s(n)  = 0 for  n<0  and  n>N.  These  limits 
allow  c(i,j)  to  be  simplified  as 

oo 

c(i>  j)  = Z x(n  " i)x(n  - j)  = r(|i  - j|).  (3-11) 

n=-«o 

Thus,  although  the  error  e(n)  is  minimized  over  an  infinite  interval,  equivalent  results  are 
obtained  by  minimizing  over  [0,  N+p-1],  Using  Eq  (3-11),  the  Wiener-Hopf  equations  are 
simplified 
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Za(i)r(|i-k|)  = -r(k),  for k = l,2,...,p 

i=l 

(3-12) 

Emin  =Za(k)r(k) 

(3-13) 

k=0 


In  matrix  form,  Eq.  (3-12)  becomes 


r(0) 

r(l)  . 

• r(p- 1) 

a(l) 

'r(l)" 

r(0 

r(0)  . 

• r(p-2) 

a(2) 

= — 

r(2) 

r(p-l) 

r(p-2)  . 

. r(0)  _ 

_a(p)_ 

_r(p)_ 

The  matrix  in  Eq.  (3-14)  is  Toeplitz  and  positive  definite.  The  resulting  poles  are 
guaranteed  to  be  within  the  unit  circle.  The  autocorrelation  method  has  been  found  to 
produce  poorer  resolution  spectral  estimates  than  other  estimators.  For  this  reason,  it  is 
not  recommended  for  short  data  records  [Kay,  1988],  Levinson  proposed  an  efficient 
algorithm  to  estimate  the  model  parameters  using  the  autocorrelation  method.  In  order  to 
discuss  Levinson  algorithm,  we  need  to  define  new  terms  called  forward  and  backward 
prediction  filters.  The  previous  definition  of  linear  prediction  in  Eq.  (3-3)  corresponds  to 
the  forward  prediction.  The  backward  prediction  filter  predicts  the  current  value  using  the 
linear  combination  of  future  values.  Using  the  Gram-Schmidt  orthogonalization  of  the 
data,  the  Levinson  algorithm  estimates  the  model  parameters  by  recursively  updating  the 
parameters  from  the  first  order  predictor  to  p-th  order  predictor.  The  Levinson  algorithm 
requires  0(p2)  operations  solving  the  Wiener-Hopf  equations,  while  Gaussian  elimination 
required  0(p3)  operations.  The  Levinson's  recursive  algorithm  requires  the  ACF  values. 
The  predictor  coefficients  for  an  i-th  order  predictor  can  be  defined  in  terms  of  the 
prediction  error  of  an  (i-l)th  order  predictor  and  the  i-th  reflection  coefficient.  It  can  be 
shown  that  the  reflection  coefficient  can  be  directly  related  to  the  forward  and  backward 
prediction  errors.  The  lattice  method  is  a variation  of  the  Levinson  algorithm.  In  the 
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lattice  method,  the  reflection  coefficient  is  calculated  in  terms  of  the  forward  and 
backward  prediction  errors.  The  lattice  filter  has  the  useful  property  that  the  coefficients 
of  the  filter,  that  is  reflection  coefficients,  are  all  less  than  one.  This  property  becomes 
important  when  the  filter  coefficients  need  to  be  quantized  for  storage  or  transmission. 
The  covariance  method  [Markel  and  Gray,  1976]  is  defined  by  setting  a = p and  b = 

N-l, 


N-l 

c(i,j)  = 22x(n-i)x(n- j)  (3-15) 

n=p 

so  that  the  error  is  minimized  only  over  the  interval  [p,  N-l],  and  all  N data  are  used  in 
calculating  the  covariance  matrix  elements  c(i  j).  In  matrix  form,  Eq.  (3-9)  becomes 


'c(l,l) 

c(2,l) 

c(l,2)  . 
c(2,2)  . 

O O 

/"“N  /-"N 

-- 
i—  *T3 

'w'  'w' 

1 

a(l) 

a(2) 

c(1.0) 

c(2,0) 

_c(p,l) 

c(p,2)  . 

• c(p>p)_ 

.a(P). 

_c(p»o) 

The  matrix  in  Eq.  (3-16)  is  symmetric  and  positive  semidefinite.  The  estimated  poles  are 
not  guaranteed  to  lie  within  the  unit  circle.  The  covariance  method  can  be  derived  from 
the  formulation  of  Prony's  method,  which  will  be  discussed  later.  For  pure  sinusoidal  data, 
the  covariance  method  yields  the  exact  locations  of  the  poles  if  we  know  the  number  of 
poles.  For  large  data  records  in  which  N » p,  the  performance  of  the  covariance  method 
is  similar  to  the  autocorrelation  method. 

The  modified  covariance  method  estimates  the  predictor  coefficients  by  minimizing 
the  average  of  the  estimated  forward  and  backward  prediction  errors.  This  method  is  also 
called  forward-backward  or  least  squares  methods  [Ulrych  and  Clayton,  1976],  The 
modified  covariance  method  is  statistically  stable  compared  to  other  methods  [Kay,  1983] 
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3.3.3.  Other  Formulations 

There  is  essentially  a limitless  possibility  for  formulations  with  various  assumptions 
and  motivations  that  finally  result  in  the  similar  implementation.  In  the  statistics  literature, 
the  term  AR  parameter  estimation  implies  the  method  to  estimate  model  parameters  by  the 
Yule-Wallcer  equations.  For  wide  sense  stationary  time  series,  the  Yule-Walker  equations 
define  the  relationship  between  the  model  parameters  and  the  ACF  of  the  time-series.  If 
we  write  the  Yule-Walker  equations  in  a matrix  form. 


T(0,0)  T(0,1)  ...  T(0,p)' 

' 1 ^ 

V' 

T(1.0)  T(l.l)  ...  T(l.p) 

a(l) 

= 

0 

,T(p,0)  T(p,l)...T(p,p), 

r 

w . 
-o'  : 

where  T denotes  ACF  of  the  time  series  and  (a(l), ...  , a(p),  a2}  are  the  AR  model 
parameters.  However,  when  limited  observation  is  available,  we  have  to  estimate  the 
ACF.  There  are  four  different  ways  to  estimate  the  ACF  from  the  time-series. 

The  autocorrelation  method  makes  the  assumption  that  data  outside  [0,  N-l]  are  zero, 
which  implies  a window  function.  The  resulting  T matrix  for  the  autocorrelation  method 
can  be  written  by 
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' T(0, 0)  T(0, 1)  ...  T(0,p)  ' 

T(1,0)  T(l»l)  •••  T(l»p) 


V T(p,0)  T(p,l) 

••  T(p,p), 

' x(0) 

0 

0 

rx(0)  x(l)  ... 

x(N-l)  0 

0 ' 

x(l) 

x(0) 

0 

0 

x(0) 

0 x(0)  ... 

x(N-l) 

0 

x(N-l) 

^ 0 ...  0 

x(0) 

x(N-l), 

0 

x(N-l) 

, o 

0 

0 x(N-l) 

(3-18) 

The  autocorrelation  method  is  also  called  the  Yule-Walker  method.  This  is  due  to  the 
equivalence  of  the  conventional  autocorrelation  estimation  for  the  Yule-Walker  equations. 
For  the  covariance  method,  the  ACF  is  estimated  only  with  the  observed  data. 


rT(0,0)  T(0,1)  ...  T(0,p)  > 
T(1,0)  T(l,l)  ...  T(l,p) 

,T(p,0)  Tfe,!)...^), 


' x(p) 

x(p-l) 

x(p  + l)  . 
x(p)  . 

x(N-l)  > 
x(N-2) 

' x(p) 
x(p  + l) 

x(p-l) 
x(p  + 2)  . 

x(0)  ' 

x(l) 

< x(0) 

x(l) 

• x(N-p-l)  , 

,x(N-l) 

x(N  - 2)  . 

. x(N-p-l)  J 

(3-19) 

We  can  expand  this  approach  to  estimate  the  ACF  assuming  zero  before  or  after  the 
observed  data.  These  methods  are  called  prewindowing  and  postwindowing,  respectively 
[Tsui,  1989], 
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The  maximum  likelihood  method  assumes  that  the  signal  is  a stationary  Gaussian 
process.  The  main  idea  of  the  maximum  likelihood  method  is  that  if  certain  values  of  the 
random  variables  are  observed,  the  parameters  which  produce  the  observation  most  likely 
are  those  which  maximize  the  probability  density.  It  can  be  shown  that  the  approximate 
maximum  likelihood  estimator  of  AR  parameters  is  found  by  solving  the  Yule-Walker 
equations  using  a suitable  estimate  of  the  ACF  [Jenkins  and  Watts,  1968],  Two  different 
ways  to  estimate  the  ACF  result  in  the  autocorrelation  and  covariance  methods. 

In  the  Maximum  Entropy  spectral  estimation,  an  explicit  extrapolation  of  a segment 
of  a known  ACF  is  performed  for  the  unknown  ACF  [Burg,  1975],  With  this  known 
ACF,  the  Maximum  Entropy  spectral  estimation  is  equivalent  to  the  AR  spectral 
estimator.  This  equivalence  is  maintained  only  for  a Gaussian  process  and  known  ACF 
values.  However,  the  data  obtained  from  most  experiments  are  a series  of  real  or  complex 
values  as  a function  of  time.  In  other  words,  we  know  only  the  data  of  the  time  series  but 
not  the  autocorrelation  of  the  input  signal.  The  ACF  calculated  using  the  data  in  the  time 
series  are  not  the  true  values  but  only  estimates.  Thus,  the  term  maximum  entropy  method 
is  no  longer  popular,  instead,  the  Burg  method  is  used.  The  Burg  method  is  quite  similar 
to  the  modified  covariance  method.  In  both,  the  sums  of  the  squares  of  the  average  linear 
prediction  error  are  minimized.  The  only  difference  between  the  two  methods  is  that  the 
Burg  method  uses  a constraint  in  the  process  of  minimization  to  assure  that  the  filter  is 
stable.  The  Burg  method  uses  recursive  procedure  to  find  AR  parameters,  which  can  be 
implemented  efficiently.  High  spectrum  resolution  with  short  analysis  interval  can  be 
achieved  with  the  method.  However,  it  has  been  known  that  spectral  estimation  depends 
on  the  initial  phase  of  the  input  signal  and  the  length  of  the  data  [Chen  and  Stegen,  1974], 
Modified  covariance  method  is  least  dependent  on  phase.  Another  problem  of  the  Burg 
method  is  spectral  line  splitting.  The  spectral  line  splitting  means  that  spectral  estimation 
shows  two  very  close  spectral  peaks  when  only  one  peak  is  present  [Fougere  et  al.,  1976], 
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In  the  correlation  matching  formulation  [Markel  and  Gray,  1976],  a match  between 
the  ACF  of  the  input  sequence  and  the  unit  sample  response  of  an  AR  model  is  desired  at 
as  many  points  as  possible.  The  estimate  of  the  ACF  matches  the  known  ACF  up  to  lag  p 
and  the  remaining  samples  are  extrapolated  such  that 


where  the  hat  denotes  the  estimated  value.  The  constraints  in  Eq.  (3-20)  leads  to  the 
autocorrelation  method. 

Prony's  method  assumes  that  the  observed  signal  is  a linear  combination  of  P complex 
exponential.  In  the  framework  of  system  theory,  this  signal  can  be  interpreted  as  an 
impulse  response  of  a linear  time-invariant  system.  The  z-transform  of  the  impulse 
response  is  called  the  transfer  function.  The  transfer  function  may  have  zeros 
incorporating  the  effects  of  initial  conditions.  Prony's  method  attempts  to  find  the  poles  of 
the  transfer  function.  In  Prony's  method,  the  process  of  finding  the  poles  of  the  system  is 
precisely  the  same  as  the  covariance  method.  Therefore,  the  covariance  method  can  have 
a new  interpretation  that  it  attempts  to  approximate  a given  data  by  a linear  combination 
of  complex  exponentials. 

A spectral  flatness  measure  (SFM)  is  a quantitative  measure  of  whiteness  (or  flatness) 
of  a spectrum.  It  can  be  shown  that  minimizing  the  prediction  error  is  equivalent  to 
maximizing  the  SFM  of  prediction  error  sequence  [Gray  and  Markel,  1974],  In  their  study 
of  the  applications  of  the  SFM  to  speech  signal  processing.  Gray  and  Markel  [1974]  used 
the  SFM  defined  as 


TxxOO, 

f"(lH  -±a(l)f„(k-l)J 


for  0 < k < p 
for  k > p. 


(3-20) 


1=1 
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S(X)-exp[}v(e)|^]  (3-21) 

-71 

where  V(0)  = ln(|x(eJ0)|2/rxx(O))  , is  the  normalized  log  spectrum  of  the  time  series  x(n). 

With  this  normalization  the  SFM,  S(X),  will  lie  between  zero  and  one,  and  equal  one  only 
for  a perfectly  flat  spectrum.  The  SFM  gives  more  weights  to  the  positive  spectral 
deviation  than  the  negative  deviation  from  a flat  average  spectrum.  In  addition  to  this, 

SFM  has  very  useful  properties  which  can  be  applied  the  AR  spectral  estimation  of  non- 
AR  processes.  There  is  a relationship  between  the  input  and  output  spectral  flatness  of  an 
AR  filter,  i.e., 

E(E)  = E(X)-^^-  (3-22) 

rJO) 

where  E and  X are  the  PSD  of  the  prediction  error  and  the  observed  data,  and  r^O),  and 
rK(0)  represent  the  powers  of  X and  E,  respectively.  The  autocorrelation  method  attempts 
to  minimize  the  error  power,  which  is  equivalent  to  minimize  the  spectral  flatness  of  the 
prediction  error.  With  the  SFM,  they  investigated  several  design  considerations  of  the 
autocorrelation  method  of  linear  prediction  of  speech:  window  types,  window  length, 
sampling  frequency  and  preemphasis  of  speech  signal.  Hi-conditioning  of  the 
autocorrelation  matrix  is  directly  related  to  the  SFM.  In  order  to  numerically  evaluate  the 
SFM  of  the  prediction  error  of  a speech  signal,  they  suggested  the  discrete  Fourier 
transform  of  twice  as  many  points  as  the  data  sequence. 

3,4.  Model  Order  Selection 

It  is  important  to  determine  the  order  of  an  AR  model.  However,  when  we  don't  have 
detailed  information  of  the  input  signal,  we  must  guess  the  order  of  the  model.  An 
intuitive  approach  to  choosing  the  model  order  is  to  monitor  the  prediction  error  as  we 
increase  the  model  order.  However,  the  minimum  prediction  error  usually  decreases 
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monotonically  with  increasing  model  order,  and  there  is  no  minimum.  If  the  order  chosen 
is  too  low,  the  spectrum  calculated  will  be  smooth.  A model  with  too  high  an  order  will 
generate  spurious  peaks,  which  may  be  confused  as  true  signals.  If  the  number  of 
sinusoidal  signal  is  known,  the  order  can  be  chosen  accordingly.  In  this  section,  some 
criteria  to  determine  the  AR  model  order  will  be  introduced. 

Akaike  suggested  two  approaches  to  choose  the  AR  model  order  [Akaike,  1970],  In 
the  first  method,  the  final  prediction  error  (FPE),  we  find  the  optimal  model  order  by 
minimizing  the  FPE,  which  is  defined  as 


FPE(p)  = 


gp(N  + p + l) 
N-p-1 


(3-23) 


where  c2p  is  the  error  power  for  p-th  order  and  N is  the  total  number  of  data  points. 

In  the  second  method  suggested  by  Akaike,  the  Akaike  information  criterion  (AIC), 
we  find  p by  minimizing  the  AIC,  which  is  defined  as 


AIC(p)  = ln(a;)+|p.  (3-24) 

Another  approach  proposed  by  Parzen  [1974],  the  criterion  autoregression  transfer 
(CAT),  is  quite  similar  to  the  two  methods  just  mentioned.  To  obtain  the  optimal  model 
order  we  minimize  CAT,  which  is  given  by 


CAT(p)  = — 
N 


z 


N-i 

No? 


N-p 
NoJ  ‘ 


(3-25) 


Other  methods  also  may  be  used  to  choose  the  optimum  order  of  the  AR  model.  For 
example,  Ulrych  and  Clayton  [1976]  suggest  that  when  the  model  order  is  in  the  range  of 
N N 

— < p < — , satisfactory  results  can  be  obtained. 

J £ 
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3.5.  Other  Factors  to  he  Considered 

There  are  several  methods  currently  used  to  estimate  the  AR  power  spectral  density 
(PSD).  They  include  the  autocorrelation,  variance,  modified  covariance,  and  Burg 
method.  Their  performance  for  large  data  records  is  comparable,  However,  for  short  data 
records  some  marked  differences  in  their  performance  exist  between  the  various  methods 
[Kay,  1988],  The  factors  that  affects  the  performance  of  the  spectral  estimators  are 

- Sampling  rates 

- Choice  of  analysis  interval 

- Window  type 

- Preprocessing. 

The  PW  Doppler  signal  of  interest  is  an  inherent  discrete  time  series  with  the  sampling 
rate  of  10.5  KHz.  However,  we  cannot  access  this  digital  signal  directly  but  only  an 
analog  signal  obtained  from  the  digital  signal.  In  our  study,  the  Doppler  signal  is  recorded 
with  the  sampling  rate  of  1 1 KHz.  The  choice  of  the  analysis  interval  is  determined  based 
on  a compromise  between  the  desire  to  have  a stable  spectral  estimate  while  minimizing 
the  averaging  of  the  time- varying  signal.  Market  and  Gray  [1976]  suggested  that  for 
speech  processing  windowing  in  both  autocorrelation  and  covariance  methods  is  advisable 
when  they  are  used  globally,  but  windowing  in  the  covariance  method  used  locally  is 
inadvisable.  Even  though  their  research  application  is  on  speech,  the  suggestion  will  be 
applicable  to  general  application  of  ARSE.  In  general,  AR  spectral  estimation  is  sensitive 
to  noise.  If  the  signal-to-noise  ratio  (SNR)  is  low,  the  result  of  the  AR  spectral  estimation 
is  no  better  than  that  of  periodogram.  At  low  SNR,  the  AR  model  is  no  longer  an  all  pole 
model.  Kay  [1988]  suggested  four  approaches  to  overcome  the  problem. 

1)  prefiltering  the  data  to  reduce  the  observation  noise. 

2)  use  ARMA  model 

3)  compensating  the  AR  parameters  for  the  biasing  effect  of  noise 

4)  use  a large  AR  model  order. 


CHAPTER  4 


PARAMETRIC  MODEL  OF  THE  DOPPLER  SIGNAL 
4.1.  Introduction 

As  stated  in  the  previous,  in  order  to  take  advantage  of  the  parametric  spectral 
estimation,  we  must  begin  with  an  accurate  structure  of  a model  that  represents  the  signal 
generating  process.  In  this  chapter,  a parametric  representation  of  the  pulsed-wave  (PW) 
Doppler  signal  is  developed.  The  block  diagram  of  the  PW  Doppler  system  is  shown  in 
Fig.  1-1.  The  objective  of  modeling  is  to  mathematically  represent  each  subsystem  and  the 
signals  throughout  the  subsystems  under  appropriate  assumptions. 

A list  of  assumptions  is  presented  in  the  next  section.  The  implications  and  limitations 
of  these  assumptions  follow.  These  basic  assumptions  will  allow  us  to  formulate  the 
Doppler  signal  for  a single  scatterer.  Only  the  main  results  of  a single  scatterer 
formulation  will  be  briefly  discussed.  Details  of  the  formulation  are  mentioned  in 
Appendix  A.  Next,  the  sample  volume  is  divided  into  subvolumes.  The  necessary 
assumptions  and  its  implication  and  limitations  are  discussed.  Subsequently,  a parametric 
model  of  a Doppler  signal  is  developed.  The  input  of  the  model  is  a white  Gaussian  noise 
and  its  output  is  pulse  wave  (PW)  Doppler  signal.  Since  the  scattering  of  ultrasound  by 
blood  is  a random  phenomenon,  all  the  information  of  the  Doppler  signal  is  specified  in 
terms  of  statistical  quantities.  In  the  last  section,  the  statistics  of  the  Doppler  signal  are 
derived. 
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4.2.  Assumptions 


Assumption  I:  Scatterers 

The  red  blood  cells  (RBCs)  are  the  main  scatterers  of  the  ultrasound.  Furthermore,  the 
movements  of  the  RBCs  are  the  same  as  the  surrounding  liquid,  called  plasma. 

Assumption  II:  Blood  as  a Random  Scattering  Medium  fRavleigh  Scattering') 

RBCs  are  identical,  independent,  point  scatterers.  Therefore,  the  scattering  is  first- 
order. 

Assumption  III:  Scattering  Process  of  an  RBC 

The  behavior  of  a red  blood  cell  is  completely  characterized  by  its  position  r(n), 
velocity  v(n),  and  complex  reflection  coefficient  A . Furthermore,  the  position  and  the 
velocity  are  discrete-time  stochastic  processes.  The  complex  reflection  coefficient  is  a 
random  variable. 

Assumption  IV:  Arterial  Blood  Flow  Conditions 

The  blood  flow  in  the  sample  volume  is  laminar,  steady  (time-invariant),  and  uniform. 

Assumption  V:  Constant  ESVDF 

The  effective  sample  volume  directivity  function  is  constant  within  a subvolume 
(ESVDF),  i.e.,  FP^ir)  = FF  = Constant,  for  each  i-th  scatterer  in  the  k-th  sample 
volume. 

Assumption  VI:  Independence  of  Amplitude  and  Phase  of  Complex  Reflection  Coefficient 
The  amplitude  and  the  phase  of  the  reflection  coefficient  are  independent. 

Furthermore,  the  phase  is  uniformly  distributed  over  [0,  27t). 

Assumption  VII:  Gaussianitv  of  k-th  Subvolume  Doppler  Signal 
ak(n)  and  3k(n)  they  are  assumed  to  be  Gaussian. 
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4.3.  Implications  and  Limitation  of  Assumptions  I-IV 

Human  blood  is  composed  of  a liquid,  called  plasma,  and  blood  cells.  More  than  99 
percent  of  the  cells  are  RBCs.  This  means  that  for  practical  purposes  the  white  blood  cells 
play  almost  no  role  in  determining  the  physical  characteristics  of  the  blood.  A human 
RBC  has  biconcave  disk  shape  with  an  average  diameter  of  8.5  urn,  a greatest  thickness  of 
2.4  urn  and  a volume  of  about  87  ± 6 nm^  [Chien,  1975],  The  average  hematocrit  of 
normal  men  is  about  42,  while  that  of  normal  women  is  about  38  [Guyton,  1986], 

Assumption  I was  experimentally  justified  by  the  fact  that  the  measured  ultrasonic 
backscattering  coefficient  from  blood  increases  almost  linearly  with  hematocrit  [Shung  et 
al.,  1984],  The  hematocrit  is  defined  as  the  percentage  of  the  blood  that  is  cells.  The 
second  part  of  assumption  I implies  that  the  measurement  of  the  velocity  of  RBCs  is  the 
same  as  the  blood  velocity. 

Assumption  II  stems  from  the  experiments  of  Shung  et  al.  [1976]  which  characterize 
the  scattering  properties  of  ultrasound  by  human  blood.  The  experiment  produced  the 
following  results 

1)  Scattering  intensity  is  proportional  to  the  hematocrit  over  the  range  from  7 to  40; 

2)  Scattering  is  isotropic,  i.e.,  the  blood  is  isotropic  scattering  medium; 

3)  Scattering  energy  increases  with  the  fourth  power  of  the  frequency. 

These  findings  are  compatible  with  the  assumption  that  the  scattering  of  ultrasound  by 
blood  obeys  the  Rayleigh  scattering  law  [Ishimaru,  1978]:  namely,  if  the  particle  size  is 
much  less  than  the  wavelength  of  the  ultrasound  (i.e.,  point  scatter)  and  if  the  scattering  is 
first  order  and  isotropic,  then  the  scattering  energy  will  increase  as  the  fourth  power  of  the 
transmitting  carrier  frequency. 

At  2 MHz,  the  wavelength  of  ultrasound  inside  the  brain  is  approximately  0.7  mm, 
which  is  about  minimum  82  times  the  length  of  the  RBC.  Therefore,  the  RBC  will  appear 
to  be  a point  scatterer.  The  assumption  of  the  point  scatter  and  the  isotropic  medium 
leads  to  the  fact  that  the  scattering  is  first-order  [Brody  and  Meindl,  1974],  The  first- 
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order  scattering  means  that  each  scatterer  is  statistically  independent  of  all  other  scatterers 
[Middleton,  1967],  In  other  words,  the  first-order  scattering  implies  that  the  scattered 
radiation  of  one  particle  does  not  interfere  with  and  is  not  rescattered  by  any  other 
particle.  Therefore,  the  Doppler  signal  from  blood  can  be  written  as  the  contribution 
(summation)  of  each  particle.  Let  D{(n)  be  the  Doppler  signal  from  a single  scatterer,  then 
the  Doppler  signal  from  the  sample  volume  D(n)  becomes 

D(n)  = X Di(">-  <4-1) 

i=sca£tcrcrs  in  a sample  volume 

The  first-order  scattering  can  also  be  supported  by  the  fact  that  the  compressibility 
and  mass  density  of  the  RBCs  differ  only  slightly  from  those  of  the  plasma.  The  scattering 
is,  therefore,  weak  and  can  be  approximated  by  first-order  scattering  [Angelsen,  1980], 

The  assumption  of  the  first  order  scattering  depends  on  the  density  of  the  scatterers. 
As  the  concentration  of  red  cells  increases,  a concentration  would  be  expected  above 
which  multiple  scattering  effects  are  prominent.  Fortunately  the  normal  density  of  the 
human  RBCs  ranges  from  38  to  42  percent  [Guyton,  1986],  and  therefore,  with  this 
justification,  the  scattering  process  can  be  assumed  to  be  first  order.  The  RBCs  tend  to 
interact  with  each  other,  e.g.  aggregation,  in  more  densely  packed  blood.  Mo  and 
Cobbold  [1986b]  proposed  a Doppler  model  for  aggregated  RBCs. 

Assumption  III  is  the  basic  abstraction  of  the  physical  properties  of  a scatterer.  This 
abstraction  is  common  to  most  mathematical  formulations  of  the  Doppler  signal.  The 
independence  assumption  now  implies  that  the  scattering  properties,  the  position,  velocity 
and  the  complex  reflection  coefficient  of  one  particle  do  not  affect  other  particles.  That  is, 
for  two  different  particles  i and  j,  we  can  write 
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P[ri(n),vi(n),Ai,rj(n),vj(n),Aj] 

= P[ri(n),vj(n),Ai]P[rj(n),vJ(n),Aj] 


(4-2) 


where  P[  ] denotes  the  probability  density  function  (PDF)  and  ri5  and  At  denotes  the 
position,  velocity,  and  scattering  coefficient  of  scatterer  i at  time  n,  respectively. 

The  reflection  coefficient  of  a single  scatterer  depends  on  the  perspective  from  the 
transducer.  The  assumption  that  the  reflection  coefficient  is  a random  variable  instead  of 
random  process  implies  that,  within  a sample  volume,  the  scatterer  does  not  change  its 
perspective  from  the  transducer.  Since  the  insonated  plane  wave  is  scattered  as  well  as 
absorbed  by  a RBC,  the  coefficient  is  expressed  in  a complex  number. 

The  identity  assumption  now  implies  that  each  scatterer  exhibits  identical  probability 
density  of  the  complex  reflection  coefficient,  that  is. 


for  all  particles  in  the  sample  volume. 

The  assumption  IV  is  the  most  ideal  condition.  However,  the  formulation  under  this 
condition  can  be  extended  to  more  complex  blood  flow.  It  is  important  to  note  that  the 
backscattering  by  blood  is  influenced  by  the  presence  of  flow  disturbance  as  well  as  the 
hematocrit.  However,  the  blood  flow  is  laminar  in  almost  all  parts  of  the  normal 
circulation  [Milnor,  1989].  The  low  velocities  and  small  diameters  encourage  laminar 
flow.  Steady  flow  implies  that  the  velocity  of  each  particle  in  a lamina  is  constant  for  all 
time.  Uniform  flow  means  the  velocity  of  each  lamina  is  the  same,  i.e.,  the  velocity  profile 
is  flat.  The  laminar  flow  condition  is  met  in  almost  all  parts  of  the  normal  circulation.  The 
only  clearly  demonstrated  exceptions  are  in  the  ascending  aorta  and  main  pulmonary  artery 
just  beyond  the  valves,  where  a brief  period  of  turbulence  may  occur  during  peak  systolic 
flow.  Pathological  deformations  of  the  arterial  structure  (e.g.,  narrowing  of  the  vascular 
lumen)  also  can  lead  to  turbulence.  [Milnor,  1989]  The  steady  flow  condition  is  not  met 
in  the  circulation  where  flow  is  pulsatile  in  a large  part  of  the  system.  Velocity  profiles 


P(A,)  = P(A) 


(4-3) 
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vary  with  time  and  are  nonparabolic  in  pulsatile  flow,  and  the  inertial  forces  of  acceleration 
are  not  taken  into  account  in  linear  equations  for  steady  flow. 

4.4.  Formulation  of  the  Doppler  Signal  by  a Single  Scatterer 
The  formulations  of  the  PW  and  CW  Doppler  signal  are  explored  by  Garbini  et  al. 
[1982],  and  Brody  and  Meindl  [1974],  respectively.  Here  we  are  interested  in  the  PW 
Doppler.  In  this  section,  we  will  briefly  discuss  the  main  results  for  the  single  scatter 
formulation.  For  more  details,  see  Appendix  A. 

If  we  express  the  transmitted  signal  as 

xT(0=ZF(t-nTip)cos(ja)0t)  (4-4) 

n=-» 

where 

F(  ):  envelope  of  a signal  pulse; 
x^:  pulse  repetition  interval; 
q0:  transmitting  angular  frequency, 

then  the  received  signal  by  a single  scatterer  at  the  location  rs  becomes 

00  7k  • r 

XR(t)  = AbG(rs)  cos(co0t-2kTs  +y)  £ F(x-nx '-)  (4-5) 

where 

Ai:  amplitude  of  the  reflection  coefficient,  \ > 0; 

y;:  phase  of  the  reflection  coefficient  0 < <j>;  < 27t; 

k = C0(/c:  transmitting  wave  number; 

c:  speed  of  ultrasound  (in  human  brain  about  1400  cm/sec); 

G(rs):  transfer  function  of  the  transducer. 

The  received  signal  is  processed  by  the  demodulator  and  subsequently  sampled.  The 
resulting  signal,  called  the  Doppler  signal,  becomes 
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XD(t)  = CAbH(rs)cos(2kTs  + Y)^8(t-nTip  - xr) 

n=-co 


where 


(4-6) 


G(r.)  F(t,  - ^ ) 

^ — (4-7) 

G(r)  F(t,  - — -)  dr 

sample  volume  0 

H(r,)  is  called  the  effective  sample  volume  directivity  function.  Using  the  Assumption  IV 
and  with  some  modification,  we  can  write  the  discrete-time  Doppler  signal  from  the 
particle,  Dt(n), 

Di(n)  = H(ri(n))  A;  cos(G>dn+Yi)  (4-8) 

where  <od  is  the  Doppler  frequency  shift. 


4,5,  Introducing  Subvolumes 

Now,  let's  divide  the  sample  volume  into  q subvolumes  as  in  Fig.  4-1.  q is  the  ratio  of 
the  total  number  of  particles  in  the  sample  volume  to  the  number  of  particles  which  leave 
the  samples  volume  during  a pulse  repetition  interval.  It  is  also  the  number  of  times  a 
given  particle  is  hit  by  a pulse.  The  number  of  q subvolumes  can  be  calculated  by 

If 

q = — ^ = xtf  (4-9) 

v 

where  v is  the  particle  velocity,  1 is  the  length  across  the  sample  volume  in  the  direction  of 
v.  fp  is  the  pulse  repetition  frequency,  and  x,  = 1/v  is  the  transit  time  of  the  particle.  For  a 
typical  2 MHz  probe,  the  v = 70  cm/sec,  1 = 7 mm  and  fp  = 10  KHz,  transit  time  is  10  msec 
and  q is  100. 
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blood  flow 


transmitted 

ultrasound 


Fig.  4-1  Dividing  sampling  volume  into  subvolumes 
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The  effective  sample  volume  directivity  function  is  defined  in  terms  of  F()  and  G(). 

F(  ) is  the  envelope  of  the  transmitted  signal  and  G(  ) represents  the  beam  pattern  of  both 
the  transmitter  and  receiver.  If  the  length  of  the  sub  volume  is  small,  then  we  can  assume 
that  the  envelop  function  is  constant  within  the  sample  volume,  that  is, 

2k  • r 

F(xr L)  » Fk , for  in  a subvolume  (4-10) 

®o 

In  addition,  it  was  found  that,  for  5 MHz  sound  wave,  the  absorption  in  blood  (hematocrit 
= 45)  is  about  0.8  dB/cm.  The  absorption  by  blood  is  approximately  linearly  proportional 
to  frequency  and  hematocrit  (for  hematocrit  < 40).  The  beam  pattern  for  a focused 
transducer  along  with  the  small  absorption  allows  us  to  assume  that,  within  a subvolume, 
each  particle  is  insonated  by  the  same  intensity  of  ultrasound,  which  justifies 

G(r)  * (F,  for  r;  in  a subvolume.  (4-11) 

therefore,  we  can  assume  that  Hk  is  a constant.  Using  the  above  assumption  of  a constant 
IF,  we  can  express  the  k-th  subvolume  Doppler  signal  at  time  n,  Dk(n),  as 


D»=Hk  £ AfWcosKn  + yfOi)) 

f=*catUnri  in  k-th  aubvolumc 

= Hkak(n)  cos(codn)+Hkpk(n)  sin(codn) 

where 


(4-12) 


a‘(n)= ^ZA?<n)  ““Crfoo) 

(4-13) 

p>)  = {XA‘(n)sin(y‘(n)). 

(4-14) 

I 


ctk(n)  and  Pk(n)  are  called  the  in-phase  and  quadrature  components  of  the  Doppler  signal 
by  the  scatterers  in  k-th  subvolume.  Note  that  in  Eq.  (4-12),  the  amplitude  and  the  phase 
of  the  reflection  coefficient  are  in  the  form  of  a stochastic  process.  This  is  because  the 
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scatterers  in  the  k-th  subvolume  are  completely  replaced  by  new  scatterers  in  a unit  time. 
Therefore,  the  reflection  coefficient  changes  with  time  so  that  now  it  becomes  a stochastic 
process.  The  complex  reflection  coefficient  of  an  RBC  depends  on  its  geometry  with 
respect  to  the  incident  plane  wave.  Since  we  assume  that  each  particle  is  independent,  the 
orientation  does  not  have  any  favorite,  which  justifies  the  uniform  distribution  of  the  phase 
over  the  interval  [0,  2n). 

As  in  Eq.  (4-13)  and  (4-14),  ak(n)  and  Pk(n)  are  sums  of  a large  number  of 
independent,  identically  distributed  random  variables.  There  are  10^  RBCs  in  a 
subvolume,  when  the  vessel  diameter  is  1 mm  and  the  length  of  the  subvolume  is  70  (xm 
and  the  hematocrit  of  40.  Therefore,  the  central  limit  theorem,  which  states  that  the  sum 
of  a large  number  of  independent  and  identically  distributed  random  variables  approaches 
to  the  Gaussian  distribution,  enables  us  to  make  the  assumption.  Since  the  Doppler  signal 
is  a linear  combination  of  the  Gaussian,  it  also  becomes  Gaussian.  Therefore  the  first  and 
second  moments  will  completely  describe  the  Doppler  signal. 

4,6,  Statistics  of  the  Doppler  Signal  from  a Subvolume 
Now,  we  will  derive  the  statistics  of  the  Doppler  signal  from  the  k-th  subvolume. 
First,  from  Eq.  (4-13)  the  mean  of  the  ak(n)  becomes 

E[ak(n)]  = ^-  Yj  E[Ak(n)  cos(yk(n))  ]. 

^ i= scatters  in  the  k-th  subvolume 

Since  each  scatterer  is  assumed  to  be  identical,  we  can  write, 

E[Ak(n)  cos(y k (n))]  = E[ A(n)  cos(y(n))]. 

Using  the  independence  assumption  of  the  amplitude  and  phase, 

a>)  =i£E[A(n)]  E[cos(Y(n))]. 

^ i 

Since  the  PDF  of  y(n)  is  uniform  over  [0  2k),  E[cos(y(n))]  = 0.  Thus,  we  have 
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Efa^n)]  = 0.  (4-15) 

We  can  easily  find  the  same  result  for  Pk(n)  such  that 

E[Pk(n)]  = 0.  (4-16) 

Now,  let's  find  the  autocorrelation  of  ak(n).  First,  consider  the  m lag,  where  m is  equal  to 
or  larger  than  one.  Since  the  particles  in  the  k-th  subvolume  is  totally  replace  by  new 
particles  and  each  particle  is  independent,  we  write 

E[ak(n)  ak(n+m)]  = E[ak(n)]  E[ak(n+m)]  = 0,  m > 0. 

When  m = 0, 

E[ak(n)  otk(n)]  = E[]T£Ak  A?  c05*^)  cos(yk)] 

> j 

A?  cos( ?*)  cos(Yk)]+E[^Ak  Ak  cos(Yk)  cos(Yf)]. 

i*j  j i 

Since  each  particle  is  independent,  the  first  term  vanishes  and  the  second  term  becomes 

N E[ A2  cos2  (y)]  = N E[ A2  ] E[cos2  (y)]  = j N 

where  a2A  is  the  variance  of  the  amplitude  and  N is  the  average  number  of  RBCs  in  the  k-th 
subvolume.  We  can  derive  the  similar  result  for  the  ACF  of  P^n).  We  can  write  for  the 
ACF  of  ak(n) 

E[  ak(n)  ak(n  + m)]  = -^Na28(m)  = oS(m),  (4-17) 

and  similarly 

E[  Pk(n)  pk(n  + m)]  = ^-Na25(m)  = a5(m).  (4-18) 

The  cross-correlation  of  ak(n)  and  Pk(n+m)  can  be  derived  in  a similar  way.  First, 
when  m is  larger  than  zero. 
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E[aKn)  P^n+m)]  = 0,  for  m>0. 

When  m = 0, 

E[ak(n)  pk(n)]  = E£SA-  Ak  cosfrf)  sin(yk)] 

» j 

= E[£Ak  Ak  cos(Yk)  sin(yk)]  + E[£Ak  Ak  cos(Yk)  sin(yf)]. 

i*J  1 

The  first  term  vanished  and  the  second  term  becomes 

NE[A’^L]=0. 

Therefore, 

Efa^n)  Pk(n+m)]  = 0,  for  all  m.  (4-19) 

4.7.  Parametric  Formulation  of  Doppler  Signal 
Since  the  flow  is  steady  and  uniform,  we  can  assume  ak(n)  = ak+1(n+l).  Therefore, 
we  can  write  ak(n)  as  a(n-k),  in  the  same  way,  Pk(n)  = P(n-k).  This  formulation  leads  to 
A(n)  and  B(n)  as  the  outputs  of  a linear  time-invariant  (LTI)  system  with  white  Gaussian 
noise  input,  i.e., 

A(n)  = b(k)  a(n  - k)  (4-20) 

lc=0 

B(n)  = Xb(k)P(n-k)  (4-21) 

k=0 

where  b(k)  = IE.  In  the  parametric  spectral  estimation  literature,  A(n)  and  B(n),  in  the 
above  form  along  with  the  white  random  input,  are  called  the  moving  average  (MA) 
model.  Even  though  this  formulation  is  implicit  in  Jones  and  Giddens,  it  will  provide 
important  insight  in  understanding  the  properties  of  Doppler  signal  from  the  viewpoint  of 
parametric  spectral  analysis.  Furthermore,  since  the  Doppler  signal,  D(n),  is  represented 
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of  the  band-limited  processes,  we  can  take  advantage  of  the  general  properties  of  the 
band-limited  processes. 

4.8.  Statistics  of  the  Doppler  Signal  from  the  Sample  Volume 
Now,  lets  consider  the  Doppler  signal  from  the  sample  volume  which  is  the 
summation  of  the  each  subvolume  Doppler  signal,  i.e., 


D(n)  = £Dk(n) 

k=0 

= 2]Hk  ak(n)  cos(codn)  + ^]Hk  pk(n)  sin(codn)  (4-22) 

k=0  k=0 

= A(n)  cos((Bdn)+B(n)sin(cDdn) 

where 

A(n)  = ]THk  ak(n),  and  B(n)  = £Hk  pk(n).  (4-23) 

k=0  k=0 

Using  the  properties  of  the  a^n)  and  Pk(n),  we  can  find  the  statistical  characteristics  of 
A(n),  B(n),  and  D(n).  First,  the  mean  of  A(n)  can  be  found  using  the  fact  that  the  a^n) 
has  zero  mean 

E[A(n)]=  E[^Hk  ctk(n)]  -£Hk  E[a‘(n)]  = 0. 

k=0  k=0 


The  autocorrelation  function  (ACF)  of  A(n)  can  be  found  by 


74 


Raa  (m) = E[A(n)  A(n  + m)] 

= ZZH'  HJE[a(n  + i)  a(n  + m + j)] 

> i 

= XZH'  HJa5(i-m  + j) 
i j 

= HjHj+mo 

where  the  value  of  Hj  is  zero  outside  the  interval  [0,  q-1].  q is  the  number  of  subvolumes. 
The  mean  and  the  ACF  of  B(n)  can  be  found  in  the  similar  way  such  that 
E[B(n)]=E[A(n)]  = 0, 

i 

Since  a^n)  and  Pk(n)  are  orthogonal  processes,  the  crosscorrelation  of  A(n)  and  B(n)  can 
be  found 


E[A(n)  B(n  + m)]  = HiHiE[a(n  + i)(3(n  + m + j)] 

> j 
= 0. 


Similarly,  we  have 


V(m)  = RAB(m)  = 0- 


ACF  of  D(n) 

R^m)  = E[D(n)  D(n+m)]  = R^m)  coscodm. 


(4-24) 


Power  spectral  Density  (PSD)  can  be  easily  found  if  we  assume  A(n)  and  B(n)  are  slowly 
varying  compared  to  cod 


SDD(co) 


Sa1(q  - op  + SAA(ca  + coH) 
2 


(4-25) 


where  S^co)  is  the  PSD  of  A(n). 


Table  4-1  Statistics  of  the  Doppler  Signal  from  the  Sample  Volume 


| Doppler  signal  representation:  D(n) 

D(n)  = A(n)  cos(coHn)  + B(n)  sin(coHn) 

Mean  of  A(n)  and  B(n) 

ErA(n)  = ErB(n)l  = 0 

ACF  of  A(n)  and  B(n) 

RAA(">)-Rm<n')  = EHkH‘**o 

k=0 

Crosscorrelation  of  A(n)  and  B(n) 

RAB(m)  = RBA(m)  = 0 

ACF  of  D(n) 

Rnn(m)  = E[D(n)  D(n+m)] 
= RA4(m)  cos(co,m) 

PSD  of  D(n) 

o (o)_SAA(©-  ©d)+SAA(o+(0d) 

2 
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4.9.  Summary 

In  this  chapter,  we  attempted  to  develop  a parametric  model  of  a PW  Doppler  signal. 
First,  all  the  necessary  assumptions  for  the  modeling  were  introduced.  By  dividing  the 
sampling  volume  into  subvolumes,  we  represented  the  scattering  characteristics  of  the 
sampling  volume  mathematically.  Equation  (4-23)  is  the  final  mathematical  representation 
of  the  Doppler  signal.  As  a result,  the  Doppler  signal  can  be  though  of  as  a frequency- 
shifted  moving-average  (MA)  process.  The  MA  process  was  determined  by  the 
characteristics  of  both  UDS  and  blood  flow.  The  parameters  of  the  MA  process  were 
determined  by  the  characteristics  of  UDS,  specifically  the  envelope  shape  of  the  insonating 
ultrasound  and  the  characteristics  of  the  transducer.  The  velocity  of  blood  determined  the 
order  of  the  MA  process.  The  variance  of  the  white  noise,  which  was  the  input  of  the  MA 
model,  was  determined  by  the  dimension  of  the  subvolume.  The  parametric  model 
developed  in  this  chapter  provides  important  insights  for  the  parametric  spectral  analysis 
of  the  Doppler  signal. 


CHAPTER  5 


DATA  COLLECTION 
5.1.  Introduction 

During  the  design  of  recording  system,  the  following  system  requirements  are 
identified.  First,  recording  Doppler  signal  for  at  least  3 hours  required  both  high  speed 
(1 1 KHz)  and  a large  storage  area  (2  bytes  * 1 1 KHz  * 3 hours  = 337.6  Mbytes).  Second, 
in  addition  to  the  Doppler  signal,  we  needed  to  record,  synchronously,  other  physiological 
data  which  were  far  slower  (about  1/100)  than  the  Doppler  signal.  For  example,  ECG 
signal  requires  86  samples/sec  with  2 byte  per  sample.  Third,  the  Streamer,  which  is  a 
commercial  software  for  data  transfer,  imposes  storage  limitation  by  the  bad  sectors  of 
hard  disk.  In  addition,  the  PDMA-16  board  generates  frequent  errors  with  the  Streamer. 

Under  these  design  specifications,  a recording  system  was  designed  with  the  hardware 
and  software  as  shown  in  Fig.  5-1. 


5.2,  Hardware 

The  recording  system  hardware  consists  of 

1)  Transcranial  Doppler  (TCD); 

2)  An  integrated  patient  monitor  (IPM); 

3)  4000A  PCM  Recording  Adapter; 

4)  H-R-D750U  Video  Cassette  Recorder; 

5)  PDMA-16  digital  interface  board; 

6)  and  IBM/PC  compatible. 
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Fig.  5-1  Recording  system  overview 
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5.2. 1.  TCP 

The  TCD  is  used  in  blood  flow  velocity  measurements  inside  the  major  arteries  of  the 
human  brain.  Table  5-1  describes  the  TCD  operating  Mode  for  blood  velocity 
measurements  in  the  middle  cerebral  artery. 

5.2.2.  IPM 

The  integrated  patient  monitor  (IPM)  provides  other  important  physiologic  signal 
such  like  electrocardiogram  (EC G),  and  invasive  arterial  blood  pressure.  We  will  refer  to 
these  signals  from  the  IPM  as  other  physiologic  signals  (OPS). 

5.2.3.  PCM  Recording  Adapter  and  VCR 

The  PCM  Recording  Adapter  contains  16  channels  for  analog  inputs  with  the  basic 
sampling  rate  of  88.2  KHz  with  word  length  of  14  bits.  Of  the  14  bits,  one  bit  is  assigned 
to  the  time  information,  and  13  bits  are  used  for  the  analog-to-digital  conversion.  The 
user  can  activate  part  of  the  16  channels  with  a predefined  sampling  frequency.  One  of  the 
main  functions  of  the  PCM  Recording  Adapter  is  to  convert  the  analog  inputs  to  13 -bit 
digital  signals.  Subsequently,  the  digital  signals  are  multiplexed  and  transformed,  or 
modulated,  into  a video  signal.  For  example,  in  case  of  2 channel  mode  where  channel  1 
and  9F  are  activated,  input  data  are  processed  in  the  sequence  of  chanl(t),  chan9F(t), 
chanl(t+l),  chan9F(t+l),  chanl(t+2),  and  so  on.  Chanl(t)  and  chan9F(t)  represent  the 
digital  data  of  channel  1 and  channel  9F  at  time  t,  respectively.  Finally,  the  VCR  receives 
this  video  signal  and  records  it  into  video  tape.  When  we  play  back  a video  tape,  two 
output  paths  are  available.  One  is  to  the  analog  output  connect.  The  other  is  to  the  digital 
output  connector.  The  digital  output  consists  of  14-bit  data  lines  and  one  data  strobe. 

The  pinout  matches  the  pinout  of  PDMA-16  interface  board  for  downloading  data  to  the 
computer.  The  PCM  Recording  Adapter  also  can  record  RS232  serial  communication 
signal.  This  signal  is  used  for  synchronization  between  recorded  data  and  important 
external  events. 
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Table  5-1  Operating  conditions  for  the  TCD 


Conditions 

Settings 

operating  mode 

Pulsed-wave  Doppler 

operating  frequency 

2 MHz 

Pulse  repetition  frequency 

10.5  KHz 

Length  of  the  sampling 
volume 

0.7  mm 
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5.2.4.  PDMA- 16  Interface  Board 

The  PDMA- 16  is  a parallel  digital  interface  board  for  data  acquisition  and  control  for 
IBM  PC/XT/AT  and  compatibles.  The  board  has  direct  memory  access  (DMA)  capability, 
which  allows  data  transfer,  either  in  8 bits  or  in  16  bits.  The  data  transfer  required  no 
intervention  of  the  central  processing  unit  (CPU). 

5.2.5.  PC/386 

The  386  PC  consists  of  80386  32-bit  CPU,  an  80387  math-coprocessor,  4 Mbytes  of 
memory,  and  a 150  Mbyte  hard  disk  with  ESDI  controller  (32  Kbyte  cache  memory). 
During  the  recording,  the  PC  provides  time  signal  through  RS232  serial  port  to  the  PCM 
Recording  Adapter.  After  the  recording,  the  data  in  a VCR  are  transfer  to  the  hard  disk  in 
the  PC.  Further  signal  processing  is  carried  out  the  PC 

5.3.  Software 

As  pointed  out  earlier,  the  commercial  software  imposed  some  difficulties.  In  order 
to  accommodate  the  recording  system  requirements,  special  purpose  software  in  assembly 
language,  was  developed  for  the  recording  system.  Figure  5-2  illustrates  the  data 
structures  used  in  the  device  driver.  Figure  5-3  illustrates  the  control  flow  of  the  program. 
The  initialization  routine  creates  necessary  data  structures  and  sets  up  the  control  registers 
in  the  direct  memory  access  (DMA)  controller.  There  are  three  main  data  structures:  a 
ring  buffer,  and  two  disk  buffers.  The  ring  buffer  holds  the  data  from  the  PCM  Recording 
Adapter  through  the  PDMA- 16  interface  board.  The  data  in  the  ring  buffer  are  decimated 
(undersampled),  and  transferred  to  one  of  the  disk  buffers.  Each  data  transfer,  are 
accomplished  by  two  processes  running  simultaneously.  When  the  VCR  is  in  play  mode, 
the  data  recorded  in  a video  tape  is  demodulated  and  transferred  to  the  PDMA- 16 
interface  board  in  the  PC.  After  the  initialization,  the  DMA  controller  transfers  the  data  to 
the  ring  buffer  as  a background  process.  The  data  transfer  to  the  ring  buffer  continues 
until  the  DMA  controller  is  reset  or  PCM  Recording  Adapter  stops  sending  data.  As  a 
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Fig.  5-2  Data  structure  for  the  device  driver 
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Initialization: 

1.  Setup  DMA  channel 

2.  Build  a ring  buffer 

3.  Build  a double  disk  buffer 


Decimation: 
Decimate  data 
and  transfer  from 
the  ring  buffer  to 
disk  buffer 
Doppler  signal:  by  4 
Other  signal:  by  1024 
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Background  process: 
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Foreground  process: 
Monitor  the  disk  buffer. 
If  the  disk  buffer  is  full, 
then  change  disk  buffer 
and  transfer  data  from 
disk  buffer  to  hard  disk 


Fig.  5-3  Control  flow  of  the  device  driver 
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foreground  process,  the  CPU  continuously  monitors  the  head  pointer  of  the  ring  buffer, 
and  transfer  new  data  to  one  of  the  disk  buffers.  During  this  data  transfer,  the  Doppler 
signal  and  other  signals  are  decimated  (undersampled)  with  preset  ratio.  For  example,  the 
Doppler  signal,  which  is  originally  sampled  with  44. 1 KHz,  is  decimated  by  the  factor  of  4. 
As  a result,  the  sampling  rate  of  the  Doppler  signal  stored  in  the  hard  disk  is  1 1025  Hz. 
Also,  the  CPU  monitors  the  amount  of  data  in  a disk  buffer.  If  the  disk  buffer  is  full,  then 
data  in  the  disk  buffer  are  transferred  to  the  hard  disk,  and  redirect  the  data  transfer  from 
the  ring  buffer  to  the  other  disk  buffer.  A BIOS  interrupt  routine,  function  2 of  interrupt 
13h,  was  used  for  the  disk  storage.  Since  the  data  are  directly  transferred  to  the  disk,  they 
are  not  retrievable  using  standard  DOS  functions.  Therefore,  special  program  was 
developed  to  retrieve  the  data  stored  in  the  disk. 


CHAPTER  6 


AUTOREGRESSIVE  SPECTRAL  ESTIMATION  OF  THE  DOPPLER  SIGNAL 

6.1.  Introduction 

In  this  chapter,  we  will  consider  the  autoregressive  spectral  estimation  (ARSE) 
methods  for  estimation  of  the  Doppler  spectrum.  The  Doppler  signal  is  collected  from  the 
Ml  section  of  the  middle  cerebral  artery.  In  chapter  3,  we  reviewed  the  general  features 
of  parametric  spectral  estimation,  specifically  autoregressive  spectral  estimation  (ARSE). 
As  the  first  step  to  apply  ARSE  to  the  Doppler  signal,  a parametric  model  for  the  Doppler 
signal  was  proposed  in  chapter  4.  However,  the  model  was  based  on  the  assumption  that 
the  flow  is  laminar,  steady,  and  uniform.  Moreover,  the  resulting  model  is  a frequency- 
shifted  MA  process  instead  of  an  MA  process.  Therefore,  this  chapter  will  begin  with 
defining  research  problems  when  we  want  to  use  ARSE  methods.  Next  we  will  present  the 
research  tools.  The  analysis  tool  will  be  used  to  find  optimal  design  factors  of  ARSE 
methods.  The  design  factors  for  the  ARSE  includes  implementation  methods,  size  of 
analysis  interval,  AR  model  order,  and  window  type. 

We  will  consider  two  implementation  algorithms  for  ARSE:  the  autocorrelation  and 
Burg  methods.  Three  analysis  intervals  of  32,  64,  and  128  points  are  used.  These 
correspond  to  2.9,  5.8  and  1 1.6  [msec].  For  these  three  intervals  and  two  algorithms,  we 
will  determine  the  best  model  order.  In  the  autocorrelation  method,  we  will  investigate 
the  effects  of  windows  such  as  rectangular,  Hamming,  and  Hanning  windows.  Table  6-1 
summarizes  the  design  considerations  and  the  analysis  procedures. 
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Table  6-1  ARSE  design  considerations  and  analysis  procedures 


Design  consideration 

Choices 

Analysis  procedures 

Algorithm 

Autocorrelation,  burg 

FPE,  AIC,  CAT 

Analysis  interval  (N) 

32,  64,  128 

Spectral  distance  measure(SDM) 

Model  order  (p) 

up  to  16 

Spectral  flatness  measure(SFM) 

For  autocorrelation 
method,  window  type: 

Hamming,  Hanning, 
Rectangular 
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6.2.  Research  Questions  for  ARSE  of  the  Doppler  Signal 
In  Chapter  4,  we  developed  a model  for  the  PW  Doppler  signal  for  steady,  laminar 
and  uniform  flow  resulting  in  a frequency-shifted  moving  average  (MA)  process. 

D(n)=  A(n)  cos(codn)  + B(n)sin(©dn)  (4-22) 

where 


A(n)-f  Hka»,  andB(n)«2Hk  P‘(n). 

k=0  k=0 


(4-23) 


The  Doppler  power  spectral  density  (PSD),  SDD(co),  can  be  found  if  we  assume  A(n)  and 
B(n)  are  slowly  varying  compared  to  ©d 


C Saa(®  - ©d>  + SaaC©  + C0d) 


(4-25) 


where  S^co)  is  the  PSD  of  A(n). 

The  number  of  q subvolumes,  or  the  order  of  the  MA  processes  in  Eq.  (4-23),  is 
calculated  by 


q= 


(4-9) 


where  v is  the  particle  velocity,  1 is  the  length  across  the  sample  volume  in  the  direction  of 
v.  fp  is  the  pulse  repetition  frequency  (PRF),  and  xt  = 1/v  is  the  transit  time  of  the  particle. 
The  ultrasound  Doppler  system  described  in  Chapter  5 is  operated  with  2 MHz  probe,  the 
length  of  sampling  volume  of  7 mm  and  pulse  repetition  frequency  of  10  KHz.  Therefore, 
from  Eq.  (4-9)  the  model  order  q is  inversely  proportional  to  the  blood  velocity.  In  most 
blood  velocity  measurements  at  the  middle  cerebral  artery,  the  maximum  peak  velocity  is 
about  150  cm/sec,  yielding  model  order  of  the  MA  process  of  47. 

Another  characteristic  of  the  model  is  that  the  coefficients  of  the  MA  process  in  Eq. 
(4-23)  are  all  positive.  Therefore,  we  can  assume  that  the  MA  process  is  low-pass.  Since 
the  Doppler  power  spectral  density  (PSD)  is  the  convolution  of  the  MA  process  and  the 
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sinusoidal  signal,  it  is  a narrowband  process  [Papoulis,  1981].  Furthermore,  as  the  model 
order  increases,  the  bandwidth  of  the  PSD  will  decrease.  If  we  use  an  ARSE  method  to 
estimate  the  PSD,  then  the  peak  frequency  of  the  spectrum  is  equivalent  to  the  location  of 
the  pole  of  the  AR  model. 

The  statement  so  far  holds  only  if  we  assume  that  the  blood  flow  is  steady,  uniform, 
and  laminar.  In  the  normal  human  artery,  blood  flow  is  laminar  with  a few  exceptions 
[Milnor,  1989].  Thus,  we  can  assume  that  the  blood  flow  at  the  middle  cerebral  artery  is 
laminar.  However,  in  most  cases,  the  Doppler  signal  is  from  pulsatile  blood  flow  is  neither 
steady  nor  uniform.  The  time-varying  nature  of  the  normal  blood  flow  implies  that  the 
Doppler  signal  is  a nonstationary  process.  However,  Mo  and  Cobbold  [1988b]  reported 
that  the  Doppler  signal  from  a normal  carotid  artery  can  be  assumed  be  wide-sense 
stationary  within  the  interval  of  10  msec.  The  diameter  of  the  normal  carotid  artery  is 
larger  than  the  middle  cerebral  artery.  Also,  the  resistance  of  vascular  bed  to  which  these 
arteries  supply  blood  is  different.  Therefore,  the  10  msec  boundary  of  wide-sense 
stationarity  may  not  hold  for  the  Doppler  signal  from  the  middle  cerebral  artery.  An 
objective  of  this  chapter  is  to  find  an  optimal  size  of  interval. 

The  nonuniform  nature  of  the  velocity  distribution  in  a sample  volume  means  that  the 
Doppler  spectral  shape  is  the  contribution  of  laminas  with  different  velocity.  This  may 
suggest  that  various  shapes  of  the  velocity  distribution  entail  ARSE  with  various  model 
orders  during  a cardiac  cycle.  Another  objective  of  this  chapter  is  to  answer  if  we  can  use 
a fixed  model  order  during  a cardiac  cycle.  In  the  uniform  blood  flow,  the  physical 
meaning  of  the  pole  directly  related  to  the  blood  velocity.  However,  this  physical  meaning 
no  longer  holds  for  the  nonuniform  velocity  distribution.  An  other  objective  of  this  chapter 
is  to  determined  physical  meaning  of  the  pole  location. 

There  are  several  ways  to  implement  an  ARSE.  These  implementation  methods  are 
discussed  in  Chapter  4.  Each  method  possess  unique  characteristics.  In  this  study,  we 
select  the  autocorrelation  method,  and  the  Burg  method.  The  main  reason  choosing  these 
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method  is  that  the  autocorrelation  is  computationally  efficient  and  suitable  for  a broadband 
process  while  the  Burg  method  is  suitable  for  a narrowband  process.  The  experimental 
results  for  both  methods  are  presented.  The  autocorrelation  method  involves  the  selection 
of  a window.  Therefore,  the  effects  of  windows  will  be  investigated. 

In  summary,  for  ARSE  of  the  Doppler  signal  from  the  middle  cerebral  artery,  we 
want  to  answer  the  following  questions 

1)  what  is  the  optimal  size  of  the  analysis  interval? 

2)  Can  we  use  a fixed  model  order?  If  we  can,  what  is  the  optimal  model  order? 

3)  What  is  the  physical  meaning  of  the  poles? 

4)  Which  method  perform  better,  the  autocorrelation  method  or  the  Burg  method? 

5)  What  are  the  effects  of  various  windows  in  autocorrelation  method? 

In  order  to  answer  the  above  questions,  we  employed  several  analysis  procedures.  The 
next  section  presents  the  details  of  these  procedures. 

6.3.  Analysis  Procedures 

6.3,1.  Best  Model  Order  Estimations 

In  order  to  investigate  the  optimal  model  order  of  the  Doppler  signal,  three  procedures  are 
used,  i.e.,  final  prediction  error  (FPE),  Akaike  information  criterion  (AIC),  and  criterion 
autoregression  transfer  (CAT).  These  procedures  are  carried  out  for  the  analysis  intervals 
of  32,  64  and  128  points.  The  details  of  these  three  procedures  are  discussed  in  Chapter 
4.  With  these  procedures,  we  want  to  answer  the  following  questions 

1)  Is  there  any  correlation  between  the  best  model  orders  and  the  cardiac  cycle?  Can 
we  use  a fixed  model  order  during  a cardiac  cycle. 

2)  Is  there  any  correlation  between  the  best  model  orders  and  the  size  of  analysis 
interval? 

3)  What  procedure  is  the  best  estimator  for  the  AR  model  order? 
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6.3,2.  Spectral  Distance  Measure  (SDNU 

A spectral  distance  measure  is  a quantitative  value  that  compares  the  similarity 
between  two  given  spectra.  There  are  many  definitions  of  spectral  distance  measure. 
Itakura  and  Saito  proposed  a spectral  distance  measure  [1968] 

d„(f(n),  f(n + 1))  = J (ev<">  - V(o  ) - 1)^  (6-1) 

where  V(w)  defines  the  difference  between  the  neighboring  two  spectra  f(n)  and  f(n+l), 

i.e., 

(6'2) 

This  measure  is  nonuniform  and  asymmetric.  Markel  and  Gray  [1976]  showed  that  the 
distance  measure  in  Eq.  (6-1)  emphasizes  the  spectral  peaks  more  than  the  valleys  or 
regions  of  low  energy  behavior.  There  are  two  reasons  for  adopting  this  measure.  The 
first  reason  is  the  nonuniform  spectral  weighting  of  the  autocorrelation  method.  That  is, 
the  autocorrelation  method  itself  weighs  the  spectral  peaks  more  heavily  that  the  valleys  in 
error  criterion.  In  fact,  the  autocorrelation  method  by  maximum  likelihood  formulation 
leads  to  the  minimization  of  the  integral  in  Eq.  (6-1).  The  other  reason  is  that,  when  the 
velocity  envelope  is  calculated  by  pole-tracing  algorithm,  we  are  more  interested  in  the 
spectral  peaks  than  the  rest  of  a spectrum.  The  physical  meaning  of  the  pole  of  an  AR 
model  is  the  dominant  lamina  of  the  blood  flow. 

Using  the  correlation  matching  properties,  we  can  rewrite  Eq.  (6-1)  as 

a2  a,TRa'  rr'2 

d(f(n),  f(n  + 1))  = — j— + In  ~ — 1 (6-3) 

a a Ra  cr  ' 

where  aT  = [1  a(l)  a(2) ...  a(p)]  is  the  impulse  response  of  the  all-zero  filter  and  R is  the 
(p+1)  x (p+1)  autocorrelation  Toeplitz  matrix.  The  quantities  associated  with  f(n)  is 
indicated  by  prime. 


91 


The  Itakura-Saito  measure  is  asymmetric.  However  we  can  define  a symmetric 
measure  by 

SDM(f(n),  f (n  + 1))  = d^f  ^ ^ + 1))  + dis(f(n  + 1)-ifCn))  (6- 

In  this  research,  the  symmetric  SDM  defined  in  Eq.  (6-4)  will  be  used.  The  main 
purpose  of  using  the  SDM  is  to  determine  an  optimal  analysis  interval.  Specifically,  we 
want  to  answer  the  following  questions 

1)  Is  there  any  correlation  between  a cardiac  cycle  and  the  degree  of  stationarity? 

2)  What  is  the  relationship  between  the  degree  of  stationarity  of  the  Doppler  signal 
and  the  size  of  analysis  interval? 

3)  Is  there  any  distortion  due  to  the  short  analysis  interval? 

6,3.3.  Spectral  Flatness  Measure  (SFM1 

The  spectral  flatness  measure  (SFM)  defined  in  Chapter  3 can  be  used  to  investigate 
the  performance  of  ARSE.  The  SFM  of  a signal  is  a measure  that  indicates  how  flat  the 
spectrum  is.  When  we  employ  an  parametric  spectral  estimation,  we  assume  that  the 
observed  data  are  an  output  of  a linear  system  excited  by  a white  noise  input.  In  linear 
prediction  formulation,  we  can  find  the  input,  or  the  prediction  error  sequence,  using  the 
inverse  filter.  By  evaluating  the  SFM  of  the  prediction  error,  we  can  determine  whether 
PSD  of  given  sets  of  data  can  be  successfully  estimated  by  an  AR  model.  Using  the 
spectral  flatness  measure,  we  want  to  answer  the  following  question. 

1)  When  we  use  a fixed  AR  model  order,  is  there  any  correlation  between  the  SFM 
and  a cardiac  cycle?  Can  we  use  a fixed  model  order  during  a cardiac  cycle?. 

2)  Is  there  any  correlation  between  the  SFM  and  the  size  of  analysis  interval? 

In  this  study,  we  will  use  a SFM  of  PSD  of  prediction  error  E(co).  For  numerical 

calculation,  we  approximate  the  integral  by  summation  using  N distinct  point  and 
rectangular  integration  with  the  result, 
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S(E)* 


exp  ^ElnlE(k)l2 

k=l 


mm1 


. k=l 


^EiE(t)r 

k=l 


~Z|E(k)|2 

iN  k=l 


(6-5) 


where  N can  be  arbitrarily  large,  and  E(k)  = E<{  exp 


j27t(k  — 1) 
N 


►,  which  is  found  by  FFT. 


In  this  form,  the  SFM  can  be  interpreted  as  a ratio  of  a geometric  to  an  arithmetic  means 
of  |E(k)|2.  The  value  of  the  spectral  flatness  measure  lies  between  zero  and  one.  A flat 
spectrum  will  yield  the  SFM  of  one. 


6,3.4.  Data  Collection 

We  use  Doppler  signal  during  two  cardiac  cycles.  The  test  Doppler  signal  is  from  the 
middle  cerebral  artery  from  a subject.  This  signal  is  collected  using  the  automatic 
recording  system  described  in  Chapter  5.  Figure  6-1  shows  the  maximum  velocity 
envelope  (MVE)  of  the  test  Doppler  signal.  The  (MVE)  is  found  by  the  percentile  method 
with  threshold  of  98  percent. 


6,4.  Results  and  Discussion 


6,4.1.  Autocorrelation  Method 

In  this  section,  we  investigated  the  performance  of  autocorrelation  method.  The 
autocorrelation  method  was  implemented  efficiently  using  the  Levinson  algorithm. 

6,4, 1, 1,  Analysis  by  best  order  estimation  procedures 
In  order  to  observe  any  correlation  between  a cardiac  cycle  and  best  model  orders, 
the  best  orders  selected  by  the  three  criteria  are  plotted  as  in  Fig.  6-2,  together  with  the 
velocity  envelope.  The  analysis  interval  has  been  chosen  as  N=64  to  compare  with  the 
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Fig.  6-1  Maximum  velocity  envelope  of  the  test  Doppler  signal 
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Fig.  6-2  Best  orders  according  to  the  three  criteria,  N=64: 
(a)  maximum  velocity  envelope,  (b)  FPE,  (c)  AIC,  (d)  CAT 
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results  provided  by  Schlindwein  and  Evans  [1990],  It  can  be  observed  that  there  is  no 
clear  relationship  between  the  best  order  chosen  by  the  three  procedures  and  the  maximum 
velocity  envelope.  This  agrees  with  the  experimental  result  by  Schlindwein  and  Evans 
[1990] 

In  order  to  investigate  the  behavior  of  the  three  procedures  according  to  various 
analysis  intervals,  we  present  histograms  of  best  order  for  three  different  analysis  intervals 
as  in  Fig.  6-3.  Values  of  mean  and  the  standard  deviation  of  the  best  order  are  given  and 
plotted  in  Table  6-2  and  Fig.  6-4. 

The  behavior  of  three  procedures  is  agreeable  up  to  the  analysis  interval  of  N=64. 
Schlindwein  and  Evans  [1990]  found  that,  the  larger  the  frame  (128  and  256  points) 
becomes,  the  better  the  performance  of  the  best  order  estimate  for  simulated  AR 
processes.  Nonetheless,  as  we  use  larger  interval  (N=128)  for  the  Doppler  signal,  the 
behavior  of  all  three  procedures  becomes  inconsistent.  This  can  be  observed  as  in  the 
wider  distribution  of  best  orders  in  the  histograms  in  Fig.  6-3,  and  the  higher  standard 
deviation  in  Fig.  6-4  (b).  One  possible  reason  for  this  inconstancy  is  the  nonstationarity  of 
the  Doppler  signal.  In  other  words,  during  a cardiac  cycle,  there  exists  various  the 
velocity  distributions  of  blood.  In  an  extreme  case,  if  we  use  a model  order  far  larger  than 
the  best  model  order,  some  spurious  peaks  can  be  observed.  The  wider  distribution  of 
best  orders  suggests  that,  for  N=128,  we  may  have  problems  using  a fixed  AR  model 
order  during  a cardiac  cycle. 

As  can  be  seen  in  Fig.  6-2,  6-3,  and  , the  behavior  of  FPE  and  CAT  are  very  similar 
independent  of  the  analysis  intervals.  However,  AIC  tends  to  estimate  higher  order  than 
the  other  two  criteria,  which  disagrees  with  the  result  by  Schlindwein  and  Evans  [1990]. 
Also,  the  standard  deviation  of  AIC  is  higher  than  those  of  other  two  criteria  independent 
of  the  analysis  intervals. 

The  means  of  the  best  order  increase  with  larger  analysis  interval  as  shown  in  Fig.  6- 
4.  This  is  also  caused  by  the  nonstationarity  of  the  Doppler  signal.  That  is,  as  we  choose 
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Fig.  6-3  Histograms  of  best  orders  (a)  FPE,  N=32;  (b)  CAT,  N=32;  (c)  AIC,  N=32; 
(d)  FPE,  N=64;  (e)  CAT,  N=64;  (f)  AIC,  N=64; 

(g)  FPE,  N=128;  (h)  CAT,  N=128;  (i)  AIC,  N=128 
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Fig.  6-4  Plots  of  (a)  means  and  (b)  standard  deviations  of  the  model  order  estimates 

according  to  analysis  intervals 
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Table  6-2  Statistics  of  the  best  model  order  estimates 
(a)  mean,  (b)  standard  deviation,  (c)  percentages  of  analysis  intervals 
less  than  and  equal  to  the  proposed  model  order. 


interval  \ criteria 

FPE 

CAT 

AIC 

32 

4.7 

4.6 

6.4 

64 

6.3 

6.2 

8.8 

128 

9.2 

9.1 

12.4 

(a) 


interval  \ criteria 

FPE 

CAT 

AIC 

32 

1.2 

1.2 

1.7 

64 

1.9 

1.8 

2.7 

128 

2.7 

2.7 

3.0 

(b) 


interval  \ criteria 

FPE 

CAT 

AIC 

32,  P<6 

86% 

86% 

45% 

64,  P<8 

89% 

90% 

54% 

128,  P<12 

95% 

95% 

60% 

(c) 
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larger  analysis  interval,  time  averaging  process  emerges.  We  could  use  smaller  interval  for 
our  analysis.  However,  a smaller  interval  may  result  in  an  unstable  estimation  of  best 
order.  In  this  study,  the  optimum  model  order  for  each  analysis  interval  is  determined  to 
be  larger  than  the  mean  of  the  estimates  from  FPE  or  CAT  by  one  or  two.  These  values 
are  P=6,  8,  and  12  for  analysis  interval  ofN=32,  64,  and  128,  respectively.  Table  6-2  (c) 
is  the  percentages  of  the  number  of  segments  that  satisfy  the  optimum  values. 

In  summary,  the  three  questions  in  section  6.2.1  are  answered.  There  is  no 
correlation  between  the  best  model  order  and  the  maximum  velocity  envelope.  For  the 
analysis  intervals  of  N=32  and  N=64,  FPE  and  CAT  performed  better  than  AIC. 

However,  when  N=128,  all  three  procedures  result  in  inconsistent  behavior.  As  we 
increase  the  size  of  analysis  interval,  the  model  order  should  increase  to  accommodate  the 
time-averaging  of  nonstationary  Doppler  signal.  When  we  chose  the  analysis  interval  of 
N=128,  we  can  no  longer  use  a fixed  model  order  during  a cardiac  cycle.  This  is  due  to 
the  larger  standard  deviation  of  three  best  orders. 

6.4, 1,2,  Analysis  by  spectral  distance  measure 

We  first  employed  6th,  8th,  and  12th  fixed  order  autocorrelation  method  for  analysis 
interval  of  32,  64,  and  128,  respectively.  These  conditions  are  determined  by  adding  1 or 
2 to  the  means  of  best  model  orders  found  in  previous  section.  Figure  6-5  shows 
continuous  Doppler  spectra  for  the  120  frames  when  N=64.  Figure  6-6  shows  the 
relationship  with  a differential  velocity  envelope  and  the  spectral  distance  measure.  For 
N=32  and  N=64,  no  clear  relationship  can  be  observed  between  the  differential  velocity 
and  the  spectral  distance  measures.  This  is  an  unexpected  result.  During  the  transition 
from  the  diastolic  mean  to  the  systolic  peak,  the  blood  is  accelerated  rapidly  as  shown  in 
Fig.  6-6  (a).  However,  the  SDM  shows  that  the  rate  of  changes  in  Doppler  PSD  is  not 
periodic.  This  may  imply  that  we  can  assume  that  the  Doppler  signal  is  stationary  for 
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Fig.  6-5  Log  scale  of  Doppler  spectra  for  N=64,  P=8:  (a)  frame  1 - 30,  (b)  frame  3 1-60, 

(c)  frame  61-90  (d)  frame  91-120 
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Fig.  6-6  Differential  maximum  velocity  envelope  (DMVE)  and  log  scale  of  spectral 
distance  measure:  (a)  DMVE  (b)  N=32,  P-6;  (c)  N=64,  P=8;(d)  N-128,  P=12 
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these  intervals.  When  N=128  as  in  Fig.  6-6  (c),  we  can  observe  some  correlation.  The 
peaks  at  0.2  and  1.1  seconds,  which  are  the  transition  from  the  diastolic  to  the  systolic, 
can  be  observed  similarity  in  the  Fig.  6-6  (c). 

We  can  observe  very  interesting  behavior  when  we  compare  the  mean  of  the  spectral 
distance  for  the  three  analysis  intervals.  The  means  of  the  spectral  distance  measures  are 
calculated  after  excluding  outsiders.  A threshold  of  95  percent  is  used  to  cut  out  the 
outsiders.  This  procedure  was  necessary  because  1)  the  values  of  a few  outsiders  are  far 
lager  than  the  other  values,  and  2)  it  is  a reasonable  assumption  that  the  Doppler  spectra 
did  not  change  drastically.  It  is  suspected  that  these  large  outsiders  are  due  to  the  ill- 
condition  of  the  autocorrelation  matrix.  The  mean  values  are  plotted  in  Fig.  6-7.  The 
mean  of  the  spectral  distance  measure  is  minimum  when  N=64.  We  may  expect  that  as  the 
analysis  interval  is  getting  the  smaller , the  mean  of  the  spectral  distance  decrease. 
However,  the  mean  when  N=32,  is  larger  than  when  N=64.  This  may  be  explained  by  the 
statistical  instability  of  the  autocorrelation  method  for  short  analysis  interval.  The  larger 
mean  value  for  N=T28,  results  from  the  nonstationarity  of  the  Doppler  signal. 

6.4, 1.3.  Analysis  bv  spectral  flatness  measure 
Figure  6-8  is  the  maximum  velocity  envelope  and  SFM  of  the  prediction  error  for  the 
three  analysis  intervals.  The  prediction  error  sequences  were  evaluated  using  the  optimal 
analysis  condition  determined  in  the  previous  section.  No  clear  relationship  can  be 
observed  between  the  velocity  envelope  and  the  SFM.  This  indicates  that  the  performance 
of  the  ARSE  is  independent  of  the  maximum  velocity  envelope  during  a cardiac  cycle. 

This  result  supports  the  possibility  of  using  a fixed  model  ARSE  for  the  Doppler  signal. 

The  mean  of  SFM  decreases  with  larger  analysis  interval  as  in  Fig.  6-9.  This  means 
that,  even  if  we  increase  the  model  order  for  larger  analysis  interval,  the  performance  of 
ARSE  deteriorates.  This  can  be  explained  by  increased  nonstationarity  of  the  Doppler 
signal.  One  more  interesting  observation  is  the  rate  of  decrease  of  the  SFM  for  larger 
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Fig.  6-7  Means  of  spectral  distance  measure 
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Fig.  6-8  Maximum  velocity  envelope  (MVE)  and  spectral  flatness  measure  (SFM):  (a) 
MVE;  (b)  N-32,  P-6;  (c)  N-64,  P-8;  (d)  N=128,  P-12 
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Fig.  6-9  Mean  of  spectral  flatness  measure  for  three  analysis  interval 
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analysis  interval.  The  performance  of  the  ARSE  for  N=32  is  far  better  than  for  N-64. 

For  N=64  and  N=128,  the  performance  of  the  ARSE  is  comparable. 

6,4, 1,4,  Effects  of  windows 

In  the  autocorrelation  method,  we  assume  that  the  signal  is  zero  outside  the  analysis 
interval.  If  no  explicit  windowing  is  carried  out,  discontinuities  at  the  beginning  and  the 
end  of  the  data  can  cause  spectral  distortion.  The  general  properties  of  windows  have 
been  covered  in  various  literature.  However,  we  interested  in  the  properties  of  windows 
when  they  apply  to  the  ARSE.  Markel  and  Gray  [1976]  demonstrated  that  a tapered 
window  can  enhance  spectral  resolution  of  a speech  signal.  We  tested  the  30th  frame  of 
the  test  Doppler  signal  with  rectangular  and  Hamming  windows  using  when  N=64  and 
p=8.  Figure  6-9  shows  the  PSD  with  the  two  different  windows.  Although  not  in  the 
figure,  the  PSD  with  Hanning  window  is  almost  identical  to  the  PSD  with  Hamming 
window.  We  can  clearly  observe  the  enhancement  of  the  spectral  peaks.  In  our  Doppler 
signal  processing,  the  spectral  peaks  represent  the  dominant  blood  flow  laminas.  There 
the  peak  enhancement  by  Hamming  window  will  allow  more  accurate  estimation  of  the 
dominant  blood  flow  laminas. 

6,4.2.  Burg  Method 

Compare  to  the  autocorrelation  method,  the  Burg  method  has  been  known  to  be  more 
stable  and  have  better  resolution  capability.  In  this  section,  similar  analysis  results  will  be 
presented  and  discussed  as  in  the  autocorrelation  method. 

6.4,2. 1.  Analysis  bv  best  order  estimation  procedures 
The  three  model  order  estimation  procedures  were  used  for  the  same  analysis 
intervals.  No  clear  relationship  was  found  between  the  best  model  orders  and  the  cardiac 
cycle.  This  agrees  with  the  experimental  result  of  the  autocorrelation  method  in  the 
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Fig.  6-10  Doppler  spectra  by  autocorrelation  method 
with  the  rectangular  and  Hamming  windows 
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previous  section.  The  histogram  is  shown  in  Fig.  6-11.  Important  statistics  are  given  in 
Fig.  6-12  and  Table  6-3. 

The  behavior  of  FPE  and  CAT  are  very  similar  independent  of  the  analysis  interval. 
However,  AIC  is  extremely  unstable.  The  degree  of  the  instability  is  higher  than  in  the 
autocorrelation  method. 

The  means  of  FPE  and  CAT  increase  for  the  larger  analysis  intervals.  The  same 
explanation  of  the  nonstationary  can  be  applied  for  this.  However,  the  standard  deviations 
of  the  three  procedures  decrease  for  larger  analysis  interval.  This  may  imply  that  the 
effects  of  shorter  interval  are  more  influential  than  the  effects  of  nonstationarity. 

There  are  some  drastic  differences  between  the  autocorrelation  and  the  Burg 
methods.  The  means  and  the  standard  deviations  of  the  Burg  method  are  larger  than  those 
of  the  autocorrelation  method.  The  larger  standard  deviations  of  the  Burg  method  suggest 
that  the  autocorrelation  method  is  more  suitable  for  the  ARSE  with  a fixed  model  order. 
Another  notable  difference  is  that  the  standard  deviations  of  the  Burg  method  decrease  for 
larger  analysis  interval,  while  in  the  autocorrelation  method  the  relationship  is  reversed. 

6,4,2.2.  Analysis  bv  spectral  distance  measure 

The  spectral  distance  measure  (SDM)  for  each  interval  is  presented  in  Fig.  6-13.  As 
in  the  autocorrelation  method,  no  clear  correlation  can  be  observed  between  the 
differential  maximum  velocity  envelope  and  the  SDM.  The  means  of  SDM  for  each 
interval  were  evaluated  excluding  the  outsiders.  The  same  95  percent  threshold  was  used 
as  in  the  autocorrelation  method.  These  means  are  plotted  in  Fig.  6-14.  The  SDM 
decreases  for  larger  analysis  interval.  In  the  autocorrelation  method,  the  relationship  is 
reversed.  A similar  interpretation  can  be  applied  to  this  as  in  the  standard  deviations  in  the 
previous  section.  That  is,  for  shorter  analysis  interval,  the  Burg  method  is  unstable 
because  of  its  short  analysis  interval. 
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Fig.  6-11  Histograms  of  best  model  order  estimates: 
(a)  FPE,  N=32;  (b)  CAT,  N=32  (c)  AIC,  N=32; 
(d)  FPE,  N=64;  (e)  CAT,  N=64;  (f)  AIC,  N=64; 
(g)  FPE,  N=128;  (h)  CAT,  N=128;  (i)  AIC,  N=128 


Model  order  Model  order 
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(a) 


(b) 

Fig.  6-12  Plots  of  (a)  means  and  (b)  standard  deviations  of  the  model  order  estimates 

according  to  analysis  intervals 


Ill 


Table  6-3  Statistics  of  the  best  model  order  estimates  (a)  mean, 
(b)  standard  deviation,  (c)  percentages  of  analysis  intervals 
less  than  equal  to  the  proposed  model  order. 


interval  \ criteria 

FPE 

CAT 

AIC 

32 

7.4 

7.2 

12.8 

64 

8.1 

7.9 

12.2 

128 

10.0 

9.9 

13.4 

interval  \ criteria 

FPE 

CAT 

AIC 

32 

3.7 

3.6 

3.6 

64 

3.0 

2.9 

3.3 

128 

2.9 

2.9 

2.5 

interval  \ criteria 

FPE 

CAT 

AIC 

32,  P<10 

79.4% 

81.1% 

35.0% 

64,  P<10 

81.1% 

83.3% 

33.6% 

128,  P<12 

79.4% 

73.2% 

16.4% 
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Fig.  6-13  Differential  maximum  velocity  envelope  and  spectral  distance  measure  (SDM)  in 
log  scale:  (a)  differential  maximum  velocity  envelope;  (b)  N=32,  P=10;  (c)  N=64,  P=10; 

(d)  N=128,  P=12 


Mean  SDM 
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Fig.  6-14  Mean  of  spectral  distance  measure  (SDM)  for  three  analysis  intervals 
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6. 4. 2. 3.  Analysis  bv  spectral  flatness  measure 
No  clear  correlation  between  the  spectral  flatness  measure  (SFM)  and  the  maximum 
velocity  envelope  as  in  Fig.  6-15.  The  means  of  the  SFM  plotted  in  Fig.  6-16.  As  in  the 
autocorrelation  method,  the  means  decrease  as  the  size  of  the  analysis  interval  becomes 
larger.  This  implies  that  the  performance  of  the  Burg  method  deteriorates  with  larger 
analysis  interval. 


6,5,  Summary 

In  this  chapter,  we  considered  the  autoregressive  spectral  estimation  (ARSE)  methods 
for  estimation  of  the  Doppler  spectrum.  The  Doppler  signal  was  collected  from  the  Ml 
section  of  the  middle  cerebral  artery.  First,  research  questions  for  applying  an  ARSE 
methods  were  identified.  The  analysis  tools  to  answer  those  questions  were  introduced 
next.  Finally,  discussion  on  the  experimental  results  were  presented  for  the  optimal 
parameters  for  ARSE.  The  following  summarize  the  answer  to  the  research  questions. 

The  choice  of  analysis  interval  was  determined  based  on  a compromise  between  the 
desire  to  have  stable  spectral  estimate  while  minimizing  averaging  of  the  time-varying 
signal.  In  the  autocorrelation  method,  the  analysis  by  the  spectral  distance  measure 
showed  that  the  analysis  interval  for  64  points  (5.8  msec)  was  the  optimum.  In  the  Burg 
method,  all  three  analyses  showed  that  the  performance  was  getting  worse  for  the  shorter 
analysis  interval.  The  estimation  instability  due  to  the  shorter  interval  overwhelmed  the 
averaging  of  the  nonstationary  Doppler  signal. 

In  both  autocorrelation  and  Burg  method,  no  clear  relationship  could  be  observed 
between  the  spectral  characteristics  and  the  cardiac  cycle.  This  results  suggested  that  a 
fixed  model  order  could  be  used  for  the  ARSE.  However,  the  autocorrelation  method 
was  more  plausible  than  the  Burg  method.  This  was  due  to  the  larger  standard  deviation 
of  the  best  model  order  estimation  procedures  in  the  Burg  method. 
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Fig.  6-15  Maximum  velocity  envelope  and  spectral  flatness  measure  of  Burg  method:  (a) 
maximum  velocity  envelope;  (b)  N=32,  P=10;  (c)  N=64,  P=10;  (d)  N=128,  P=12 


Mean  SFM 
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Fig.  6-16  Mean  of  spectral  flatness  measure  (SFM)  for  three  analysis  intervals 
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Among  the  three  best  model  order  estimation  procedures,  FPE  and  CAT  were 
comparable.  However,  AIC  exhibited  instability.  This  instability  was  more  severe  in  the 
Burg  method.  The  best  model  orders  were  determined  by  adding  1 or  2 to  the  mean  of 
FPE  or  CAT  results.  The  model  order  should  increase  for  large  analysis  interval. 

Empirical  observation  indicated  that  the  Doppler  spectrum  consisted  of  one  or  two 
spectral  peaks.  This  might  imply  that  there  were  at  most  two  dominant  laminas  of  blood 
flow.  Therefore,  we  could  conclude  that  the  physical  meaning  of  the  pole  of  an  AR  model 
was  the  dominant  lamina  of  the  blood  flow. 

The  analysis  by  the  best  model  order  estimation  procedure  indicated  that  the  Burg 
method  required  more  model  order  and  the  variance  of  the  model  order  was  larger  than 
the  autocorrelation  method.  The  analysis  by  SDM  showed  that  the  Burg  method  was  more 
unstable  for  short  data  length. 

Finally,  a taper  window  enhanced  the  spectral  resolution  of  Doppler  signal.  Therefore, 
in  order  to  find  a dominant  lamina,  using  a tapered  window  such  like  Hamming  window 
was  recommended. 


CHAPTER  7 


BLOOD  VELOCITY  ENVELOPE  ESTIMATION 
7. 1 . Introduction 

In  Chapter  6,  we  developed  optimal  design  parameters  for  the  autoregressive  spectral 
estimation  (ARSE).  In  this  chapter,  velocity  envelope  estimator,  based  on  the  ARSE,  is 
proposed.  The  algorithm  can  be  divided  into  three  main  parts  as  shown  in  Fig.  7-1.  The 
optimal  model  orders  of  the  ARSE  in  Chapter  6 are  used  for  the  velocity  envelope 
estimation.  The  pole-tracing  algorithm  traces  the  highest  pole  of  the  velocity  profile  larger 
than  a threshold  value.  The  postprocessor  is  employed  for  noise  reduction.  The 
performance  of  the  velocity  envelope  estimators  are  evaluated  by  comparing  it  with  flow 
data  obtained  by  an  electromagnetic  flowmeter.  A linear  regression  was  performed 
between  the  flow  waveform  and  various  velocity  envelope  estimators. 


7,2.  Pole-Tracing  Algorithm 

A new  method  is  proposed  in  this  study  to  find  a velocity  waveform  using  an  ARSE. 
The  maximum  velocity  envelope  traces  the  maximum  velocity  of  a Doppler  spectrum.  The 
new  method  traces  a pole  of  the  AR  model  that  satisfies  certain  conditions.  The  pole  is 
found  in  the  following  manner 

Step  1.  Find  poles  where  the  PSDs  are  larger  than  a threshold  value 

Step  2.  Among  the  poles  found  in  step  1,  select  the  pole  whose  angle  is  the  highest. 
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*ARSE:  Autoregressive  spectral  estimator 


Fig.  7-1  Overview  of  the  proposed  Doppler  signal  processing 
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A careful  examination  of  the  Doppler  spectra  reveals  at  most  two  significant  peaks. 

In  other  words,  at  most  two  pole-pairs  are  significant.  A physical  interpretation  of  these 
poles  may  be  that  the  blood  flow  consists  of  two  significant  laminas  of  constant  velocities 
and  the  velocities  are  interpreted  as  the  poles.  In  our  method  these  poles  are  selected  by 
comparing  the  PSD  at  the  poles  with  a threshold  value.  The  evaluation  of  the  PSD  is 
carried  out  without  the  gain  factor.  That  is,  the  PSD  of  the  Doppler  spectra  excited  by  a 
white  noise  with  a unit  variance.  This  feature  allows  the  velocity  envelope  estimator  to  be 
adaptive  to  the  various  power  levels  during  a cardiac  cycle.  In  our  method,  the  threshold 
value  was  set  to  10  dB.  These  poles  are  then  sorted  and  the  highest  pole  is  selected.  If 
there  are  no  poles  above  the  threshold  value,  then  the  5 dB  is  subtracted  from  the  previous 
threshold  level  until  a pole  is  found. 

Figure  7-2  is  the  pole-tracing  waveform  (PTW)  based  on  the  autocorrelation  method 
and  Burg  method  with  the  optimum  analysis  conditions.  For  N=32  and  N=64,  the  PTW  is 
very  noisy.  This  is  due  to  the  instability  of  the  autocorrelation  method  for  N=32.  This 
instability  was  observed  by  the  analysis  of  spectral  distance  measure.  As  the  size  of 
analysis  interval  grows,  the  waveform  becomes  smoother. 

7,3,  Postprocessor 

When  we  use  smaller  analysis  interval,  the  PTW  becomes  quite  noisier.  This  noise  is 
generated  by  the  inherent  nature  of  the  ARSE.  When  we  employed  pole-tracing  algorithm 
for  short  analysis  interval  (N=32,  and  N=64),  we  observed  noisy  waveform.  A post 
processor  is  proposed  to  filter  out  the  noise.  It  consists  of  a median  filter  followed  by  an 
all-zero  filter.  The  median  filter  is  a nonlinear  filter  and  its  properties  are  studied  by  many 
researchers  [Beaton  and  Tukey,  1974],  The  median  filter  has  been  used  for  the  smoothing 
of  data  with  the  following  characteristics  [Rabiner  and  Schafer,  1978] 

1)  The  data  have  sharp  discontinuities  of  reasonable  duration. 
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Fig.  7-2  Maximum  velocity  envelope  (MVE)  and  pole-tracing  waveforms  (PTW)  based  on 
autocorrelation  method:  (a)  MVE;  (b)  PTW,  N=32,  P-6;  (c)  PTW,  N-64,  P=8;  (d)  PTW, 

N=128,  P-12 
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2)  The  data  have  roughness  due  to  noise  in  the  measurement  (the  noise  may  be 
inherent  in  the  measurement  process  itself  being  imperfect) 

3)  The  data  have  sharp,  isolated  discontinuities  of  very  short  duration  which  are  to  be 
eliminated.  Such  discontinuities  may  be  due  to  imperfect  analysis  procedures. 

These  characteristics  of  the  median  filter  are  ideal  when  we  process  a signal  from  human 
body  such  as  speech  and  blood  velocity.  We  used  a running  3 point  median  and  a 3 point 
Hanning  filter  with  coefficients  of  0.25,  0.5,  and  0.25.  Figure  7-3  shows  comparisons 
between  inputs  and  outputs  of  the  postprocessor.  As  expected  the  sharp,  isolated 
discontinuities  of  very  short  duration  are  eliminated.  Also,  the  waveform  is  smoothed 
because  of  the  lowpass  Hanning  filter. 

7.4,  Performance  Evaluation. 

In  order  to  evaluate  the  performance  of  our  method,  we  compare  the  blood  velocity 
waveform  with  the  blood  flow  waveform  measured  by  an  electromagnetic  flowmeter.  We 
measure  both  blood  velocity  and  flow  at  the  femoral  artery  of  a pig.  The  blood  flow  was 
normalized  by  its  minimum  and  maximum  value  so  that  the  value  lies  in  between  zero  and 
one.  Figure  7-4  shows  the  waveform  of  the  normalized  flow,  together  with  a velocity 
waveform.  The  blood  velocity  is  the  pole-tracing  waveform  (PTW)  based  on  the 
autocorrelation  method  with  the  analysis  conditions  of  N=64  and  P=8.  A linear  model 
was  used  to  relate  the  normalized  blood  flow  and  the  blood  velocity  such  that 

F = Po+Prv  (7-1) 

where  F and  V are  the  flow  and  velocity  vectors.  The  performance  was  determined  by  the 
variance  of  error  which  is  defined  as 
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<=S(f,-(P„+P,v,))1  (7-2) 

i=l 

where  f"  and  v,  are  the  elements  of  the  flow  and  velocity  vectors.  We  tested  the  velocity 
estimation.  We  test  the  performance  of  the  following  nine  methods 

MVE:  maximum  velocity  envelope  by  the  percentile  method  with  95  percent  of 
threshold 

PKA128:  PTW  by  autocorrelation  method  with  N=  128,  P=12 
FPKA64:  postprocessed  PTW  by  autocorrelation  method  with  N=64,  P=8 
FPKA32:  postprocessed  PTW  by  autocorrelation  method  with  N=32,  P=6 
PKB128:  PTW  by  Burg  method  withN=128,  P=12 
FPKB64:  postprocessed  PTW  by  Burg  method  with  N=64,  P=10 
FPKB32:  postprocessed  PTW  by  Burg  method  withN=32,  P=10 
CA128:  mean  velocity  by  complex  autocorrelation  method  N=128 
CA64:  mean  velocity  by  complex  autocorrelation  method  N=64 
The  details  of  the  complex  autocorrelation  were  presented  in  Chapter  2.  Figure  7-5  is  the 
flow  and  the  mean  velocity  by  the  complex  autocorrelation  method.  The  variances  of 
error  for  various  methods  are  given  in  Fig.  7-6.  The  postprocessed  PTW  with  N=64  is  the 
best  method.  The  PTW  by  autocorrelation  method  is  better  than  the  PTW  by  Burg 
method. 


7,5,  Summary 

In  this  chapter,  a blood  velocity  envelope  estimation  based  on  the  autoregressive  spectral 
analysis  methods  was  described.  The  physical  interpretation  of  the  velocity  waveform  is 
velocity  of  the  fastest  moving  laminar  assuming  that  the  velocity  distribution  consists  of  at 
most  two  laminas  of  constant  velocity.  The  autocorrelation  method  and  the  Burg  method 
were  used.  The  design  factors  in  Chapter  6 such  as  the  model  order  and  the  size  of 
analysis  interval  were  used  with  the  various  techniques. 
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Data  points 

Fig.  7-4  Normalized  flow  and  velocity  after  linear  regression:  The  velocity  is 
postprocessed  pole-tracing  waveform  by  autocorrelation,  N=64,  P=8 
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Fig.  7-5  Normalized  flow  and  velocity  after  linear  regression:  The  velocity  is  found  by 

complex  autocorrelation  method 


127 


Analysis  methods 


Fig.  7-6  Variance  of  error  for  various  velocity  estimation  techniques. 
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The  performance  of  each  method  was  determined  by  the  comparison  with  the  flow 
data  measured  by  an  electromagnetic  flowmeter.  In  case  of  autocorrelation  method,  best 
performance  was  for  N=64  and  P=8.  This  result  agreed  with  the  analyses  in  Chapter  6. 
The  Burg  method  was  found  to  be  inferior  to  the  autocorrelation  method.  We  also 
evaluated  the  performance  of  a time-domain  method,  namely  the  complex  autocorrelation 
method.  Even  though  the  complex  autocorrelation  method  was  an  efficient  method,  the 
performance  was  at  least  comparable  with  the  maximum  velocity  envelope  estimation 
method. 


CHAPTER  8 


CONCLUSION 

The  objective  of  this  research  was  twofold.  One  was  to  develop  a parametric  model 
of  the  pulse-wave  ultrasound  Doppler  signal.  The  other  was  to  propose  a blood  velocity 
estimation  technique  based  on  the  autoregressive  spectral  estimation  (ARSE).  Under  ideal 
assumptions,  the  Doppler  signal  was  modeled  as  the  output  of  a frequency-shifted  moving 
average  (MA)  process.  The  variance  of  the  inputs  and  the  system  coefficients  were 
determined  by  the  characteristics  of  both  ultrasound  Doppler  system  and  the  blood  flow. 
This  could  provide  the  insight  in  designing  a autoregressive  spectral  analysis  of  the 
Doppler  signal.  In  order  to  analyze  the  Doppler  signal,  data  acquisition  system  was 
designed.  Two  experiments  were  conducted  for  the  Doppler  signal  from  the  middle 
cerebral  artery  of  a normal  subject  and  from  the  femoral  artery  of  a pig.  The  first 
experiment  considered  the  design  problem  when  we  want  to  apply  ARSE  methods  to  the 
Doppler  signal.  In  the  second  experiment,  we  proposed  a blood  velocity  envelope 
estimator  and  evaluate  the  performance. 

The  first  experiment  began  with  identifying  the  problems  when  we  want  to  apply 
autoregressive  spectral  estimation  (ARSE)  methods.  Five  analysis  tools  were  employed  to 
solve  these  problems:  1)  final  prediction  error  (FPE),  2)  criterion  autoregression  transfer 
(CAT),  3)Akaike  information  criterion  (AIC),  4)  spectral  distance  measure  (SDM),  and  5) 
spectral  flatness  measure  (SFM).  Two  algorithms  (autocorrelation  and  Burg  methods) 
were  tested.  Three  sizes  of  analysis  interval  (N)  were  test:  32,  64  and  128  points.  These 
correspond  to  2.9,  5.8,  and  11.6  [msec]. 
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The  choice  of  analysis  interval  was  determined  based  on  a compromise  between  the 
desire  to  have  stable  spectral  estimate  while  minimizing  averaging  of  the  time-varying 
signal.  In  the  autocorrelation  method,  the  analysis  by  the  spectral  distance  measure 
showed  that  the  analysis  interval  for  64  points  (5.8  msec)  was  the  optimum.  In  the  Burg 
method,  all  three  analyses  showed  that  the  performance  was  getting  worse  for  the  shorter 
analysis  interval.  The  estimation  instability  due  to  the  shorter  interval  overwhelmed  the 
averaging  of  the  nonstationary  Doppler  signal. 

In  both  autocorrelation  and  Burg  method,  no  clear  relationship  could  be  observed 
between  the  spectral  characteristics  and  the  cardiac  cycle.  This  results  suggested  that  a 
fixed  model  order  could  be  used  for  the  ARSE.  However,  the  autocorrelation  method 
was  more  plausible  than  the  Burg  method.  This  was  due  to  the  larger  standard  deviation 
of  the  best  model  order  estimation  procedures  in  the  Burg  method. 

Among  the  three  best  model  order  estimation  procedures,  FPE  and  CAT  were 
comparable.  However,  AIC  exhibited  instability.  This  instability  was  more  severe  in  the 
Burg  method.  The  best  model  orders  were  determined  by  adding  1 or  2 to  the  mean  of 
FPE  or  CAT  results.  The  model  order  should  increase  for  large  analysis  interval. 

Empirical  observation  indicated  that  the  Doppler  spectrum  consisted  of  one  or  two 
spectral  peaks.  This  might  imply  that  there  were  at  most  two  dominant  laminas  of  blood 
flow.  Therefore,  we  could  conclude  that  the  physical  meaning  of  the  pole  of  an  AR  model 
was  the  dominant  lamina  of  the  blood  flow. 

The  analysis  by  the  best  model  order  estimation  procedure  indicated  that  the  Burg 
method  required  more  model  order  and  the  variance  of  the  model  order  was  larger  than 
the  autocorrelation  method.  The  analysis  by  SDM  showed  that  the  Burg  method  was  more 
unstable  for  short  data  length. 

Finally,  a taper  window  enhanced  the  spectral  resolution  of  Doppler  signal.  Therefore, 
in  order  to  find  a dominant  lamina,  using  a tapered  window  such  like  Hamming  window 
was  recommended. 
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In  the  second  experiment,  a blood  velocity  envelope  estimation  technique  was 
proposed.  This  technique  consisted  of  the  ARSE,  pole-tracing  algorithm,  and  a 
combination  of  median  and  three-point  Hamming  filter.  The  performances  of  this 
technique  under  various  conditions  were  evaluated  by  comparing  with  the  blood  flow 
signal  measured  by  an  electromagnetic  flowmeter.  Also,  a time-domain  analysis  technique 
(complex  autocorrelation  method)  was  included  for  the  evaluation.  The  evaluation 
involved  a linear  regression  of  two  waveforms.  The  variance  of  error  was  used  for  the 
quantitative  evaluation.  As  a result,  the  best  condition  was  found  to  be  the 
autocorrelation  method  with  N=64  and  P=8.  This  conforms  the  experimental  result  in  the 
first  experiment. 

In  summary,  three  major  accomplishments  were  achieved  in  this  study,  First,  we 
showed  that  the  PW  Doppler  signal  could  be  modeled  as  a frequency  shifted  MA  process. 
Second,  we  determined  the  optimal  model  order  and  analysis  interval  for  autoregressive 
spectral  analysis.  Finally,  we  developed  a blood  velocity  estimation  technique  using  the 
pole-tracing  and  median  filtering.  The  autocorrelation  for  P=8  and  N=64  resulted  in  the 
best  performance. 

Some  suggestions  for  further  development  in  this  research  are  as  follows.  The  design 
criteria  found  in  this  study  can  be  applied  for  adaptive  spectral  estimation  techniques. 

Since  the  Doppler  signal  is  nonstationary  in  nature,  the  adaptive  spectral  estimation  could 
be  more  suitable  for  the  spectral  characteristics. 

One  of  the  important  parameters  that  this  study  doesn't  deal  with  is  the  variance  of  the 
white  noise  of  the  model.  By  finding  the  ratio  of  the  variance  to  the  total  signal  power,  we 
can  determine  the  length  of  the  subvolume.  The  length  of  the  subvolume  is  proportional 
to  the  blood  velocity.  One  significant  feature  of  the  ratio  is  that  it  is  the  function  of  blood 
velocity,  characteristics  of  the  transducer,  and  the  shape  of  the  transmitted  ultrasound.  A 
recent  study  showed  that  the  Doppler  signal  power  was  a function  of  many  variables 
namely,  hematocrit,  shear  rate,  and  turbulence  [Shung  et  al.,  1992],  This  hindered  the 
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mean  flow  estimation.  However,  the  ratio  is  independent  of  the  variables.  Therefore,  the 
ratio  can  be  used  to  extract  mean  flow  information. 


APPENDIX 


DOPPLER  SIGNAL  FROM  A SINGLE  SCATTERER 
A.l.  Introduction 

As  stated  in  chapter  4,  we  can  assume  that  each  red  blood  cell  (RBC)  is  statistically 
independent  and  that  the  scattering  is  first-order.  These  assumptions  allow  us  to  easily 
extend  a single  scatter  formulation  to  a multiple  Scatterer  formulation.  That  is,  the 
Doppler  signal  from  multiple  scatters  is  the  linear  combination  of  the  Doppler  signal  from 
a single  scatterer.  In  this  Appendix,  a mathematical  representation  is  developed  for  the 
pulse-wave  (PW)  Doppler  signal  by  a single  point  scatterer.  Each  step  in  the  formulation 
will  provide  invaluable  insight  to  understand  the  Doppler  signal.  Brody  and  Meindl 
[1974],  Garbini  et  al.  [1982],  and  Mo  and  Cobbold  [1986b]  developed  the  similar 
representations. 

The  formulation  will  be  based  on  the  functional  diagram  in  Fig.  A-l.  The  following  is 
an  outline.  First,  a representation  of  transmitted  signal  is  presented.  The  transducer 
transforms  this  electric  signal  into  an  acoustic  wave.  A particle  scatters  this  wave  and  the 
expression  for  the  backscattered  wave  is  derived.  The  same  transducer  receives  the 
backscattered  wave  and  transforms  it  into  an  electric  signal.  The  received  signal  is 
processed  by  a mixer  and  low-pass  filter  (demodulator)  and  a sampler.  Finally,  an 
expression  of  the  Doppler  signal  is  derived  for  a single  scatterer. 
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carrier 
exp(]  cjot) 


Fig.  A-l  Functional  diagram  of  TCD 
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A.2.  Transmitted  Signal 
Let  a single  transmitted  pulse,  X,.0(t),  be  of  the  form 

XT0  (t)  = Re[XTo(t)]  = Re[F(t)exp(jco0t)]  (A-l)  • 

where  F(t)  represents  the  envelop  and  exp(jco0t)  represents  the  carrier.  Throughout 
Appendix  A,  the  topbar  indicates  a complex  variable.  Given  the  expression  for  a single 
pulse,  we  can  represent  the  transmitted  pulse  train,  XT(t),  as  the  convolution  of  the  single 
pulse  with  a train  of  unit  impulses  with  a pulse  repetition  period  of  x^See  Fig.  A-2) 


Xx(t)= 

-oo  n=-oo 

= i;xT.(t-nT,). 

n=-oo 

If  we  assume  that 


27tm 


where  m is  an  integer,  then  the  transmitted  pulse  train  can  be  written  as 


Xx(t)  = ^F(t-nTq))exp(jco0t). 

n=oo 


(A-2) 


A.3.  Transmitted  Wave 

We  assume  that  the  transducer  transforms  an  electric  signal  into  an  acoustic  plane 
wave  with  the  transfer  function  denoted  by  GT.  That  is,  the  acoustic  plane  wave 
propagating  in  the  i direction  at  the  position  r is  given  by 


voltage 
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Fig.  A-2  Transmitted  signal:  pulse  repetition  frequency  (PRF)=10.41  KHz 
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PT(t,r)  = GT(i)XT(t-— ) (A-3) 

©0 

where 

A* 

i : unit  vector  in  the  direction  of  incident  plane  wave 
k = cOg/c:  transmitting  wave  number 

A A. 

GT(i  ):  transmitter  transfer  function  in  the  direction  i 
c:  speed  of  sound  (in  human  brain  about  1400  m/sec). 

The  bold  characters  represent  vectors,  and  the  hat  represents  a unit  vector.  Figure  A-3 
illustrates  the  coordinate  system  and  important  vectors  throughout  the  analysis. 


A.4.  Scattered  Wave 

We  consider  a single  scatterer  to  be  located  at  rt.  When  we  observe  the  scattered 
wave  in  a far  distance,  the  scattered  wave  behaves  as  a spherical  wave.  We  define  r1  as  a 
new  coordinate  whose  reference  point  is  in  the  center  of  the  scatterer.  The  scattered 
wave,  P5(t,r'),  is  given  by 


r con 


(A-4) 


where  A (6 , i ) is  called  the  complex  scattering  coefficient.  The  complex  scattering 
coefficient  represents  the  amplitude  and  phase  of  the  scattered  wave  in  the  direction  of  6 
when  the  particle  is  insonated  by  a plane  wave  propagating  in  the  direction  of  i with  unit 
amplitude.  Replacing  r'  using  the  relation  between  r1  with  r,  i.e.,  r'  = r - rs , the  scattered 


wave  becomes 
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Transducer 


A 

O 


Fig.  A-3  Coordinate  system 
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P,(t,r)  = 4^PT(t-^^,0. 
k-rj  <o0 


(A-5) 


A.5.  Backscattered  Wave 

After  transmitting  a pulse,  the  transducer  waits  for  the  scattered  wave  to  return, 
which  travels  in  the  reverse  direction  of  the  incident  wave.  In  this  case,  the  returning 
wave  is  called  the  backscattered  wave.  As  a result,  the  backscattering  acoustic  wave, 
P b*(t>r)>  can  be  found  by  replacing  6 by  - i , i.e., 


lr-r.l  ©0 


(A-6) 


A.6,  Received  Signal 

At  r = 0,  the  transducer  receives  the  backscattered  wave  and  transforms  the  wave  into 
an  electric  signal,  XR(t).  Therefore,  we  can  write 


XR(t)  = GR(-i)  Pb»(t,r  = 0) 

=cR(-i)  pl(t.Mjl±Oj  r>) 


(A-7) 


A A 

where  GR(-  i ) is  the  receiver  transfer  function  in  the  - i direction.  Replacing  PT  by  Eq. 
(A-3),  the  received  signal  is  expressed  in  terms  of  the  transmitted  signal 


x„(t)=A(-i,i)  xT(t.M% 

r.  on 


(A-8) 


Let 
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G(O.GT(i)GR(-0 


(A-9) 


then  G(r,)  becomes  a function  that  characterizes  the  transducer  in  both  cases  of 
transmitting  and  receiving.  Furthermore,  we  can  write  the  complex  backscattering 
amplitude  such  that 

Abexp(jy)  = A(-i,i)  (A-10) 

where  \ is  the  backscattering  amplitude  and  y is  the  backscattering  phase.  Replacing  the 
transmitted  signal,  XT,  with  Eq.  (A-l),  and  taking  the  real  component  of  XR(t),  we  have 


XR(t)  = Re[XR(t)] 

“ ?kr 

= AbG(rs)  cos(co0t-2k-r,+y)  £ F(x-nx,p L) 

n=-oo  ® n 


(A-l  1) 


where  k = k i . 


A.7.  Mixer  Output 

The  received  signal  is  multiplied  with  a sinusoid  of  the  same  frequency  as  the  original 
transmitted  signal.  The  resulting  quantity  contains  sinusoidal  term  whose  frequencies  are 
the  sums  of  the  Doppler  shifted  frequencies.  This  signal  is  lowpass  filtered  to  remove  the 
high  frequency  terms.  The  output  of  the  filter,  XM(t),  becomes 

“ 9k . r 

XM(t)  = AbG(rJ)  cos(2k-rs  +y)  £ F(t-nx,p L)-  (A-12) 

n=-<o  © 0 

A. 8,  Doppler  Signal 

The  continuous  function  of  time  in  Eq.  (A-12)  is  sampled  at  interval  of  x We  can 
express  this  sampling  process  as  the  multiplication  of  the  mixer  output  by  a train  of 
impulse  functions.  Let  the  sampled  mixer  output  be  the  Doppler  signal,  XD(t),  then 
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xD(t)  = xM(t)  2;5(t-ttt,-dTip)  (A-i3) 

n=-<*> 

where  d is  the  range  delay  constant  0 < d < 1.  For  convenience  let  x,=  dx^.  Replacing  the 
mixer  output,  XM(t),  with  Eq.  (A- 12),  the  Doppler  signal  becomes 


" 2k  .r 

XD(t)  = Ab G(r,)  cos(2k-rI+y)^  F(t-nx 5-) 

(A-l 

x ^(t-mx.p-x,) 

m=-°o 

For  m = 1 + n and  using  the  equality  that  F(t  - a)  6(t  - P)  = F(P  - a)  6(t  - P),  the  Doppler 
signal  can  be  written  as 


XD  (t)  = Ab  G(r.)  cos(2k  • r,  + y)  £ FOx, 

1=-O0 


+ tr- 


2k-rt) 

®o 


x^6(t-nx^-xr). 

n=— oo 


(A- 15) 


If  we  consider  only  one  sample  volume  closest  to  the  transducer,  we  can  set  1 = 0 and  the 
above  equation  becomes 


2k  r 


XD(t)-Ab  G(r.)  F(xr )cos(2k  • rs  +y)  £5(t-nx  -xr). 

co 


(A- 16) 


If  we  define  H(r),  called  effective  directivity  function  such  as 


H(r.)  = G(r.)F(Tr-  ^ ) 


C 0 , 


(A- 17) 


then,  the  Doppler  signal  becomes 
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xd (0  = AbH(rs)cos(2k ■ r,  + y)  £ 8(t  - nx,,  - x, ) . (A- 1 8) 

n=-oo 

A. 9.  Summary 

Here,  a mathematical  model  of  the  Doppler  signal  by  a single  scatterer  is  developed. 

If  we  express  the  transmitted  signal  as 

00 

Xx(t)=  2F(t_nTip)  (A-l) 

n=-oo 

then,  the  received  signal  becomes 

“ 2k  r 

Xr(0  = AbG(rs)  cos(o0t-2k-rf  +y)  £F(x-nx L)-  (A-10) 

«,=-»  ©0 

The  received  signal  is  processed  by  the  demodulator  and  sampled.  The  resulting  signal, 
called  Doppler  signal,  is 

xD(0“AbH(rJcos(2k-ri  + y)£8(t-nxip  - xr).  (A- 17) 

n=*oo 

Therefore,  the  Doppler  signal  by  a single  scatterer  located  at  rs  is  expressed  by  the 
following  parameters 

1)  carrier  wave  number,  k = &Jc  where  co0  is  carrier  frequency  and  c is  the  speed  of 
sound  in  the  medium 

2)  pulse  repetition  interval,  x^ 

3)  range  delay,  xr 

4)  effective  sample  volume  directivity  function  H(r,)  which  is  the  multiple  of  transfer 
function  of  the  transducer,  G(rJ),  and  envelop  shape  of  the  transmitting  signal,  F(t) 

5)  backscattering  amplitude  and  phase  of  the  scatter,  \ and  y. 

Among  these,  the  system  parameters  of  the  ultrasound  Doppler  instrument  are  carrier 
frequency,  pulse  repetition  interval,  range  delay,  transfer  function  of  the  transducer.  The 
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single  scatterer  is  characterized  by  the  backscattering  coefficient  and  velocity  of  the 
scatterer. 
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